ex5 regression
ex5.stan

data{
  int N;
  int K;
  vector[N] y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> s;
}
model{
  vector[N] m=X*b;
  y~normal(m,s);
}
generated quantities{
  vector[N] y1;
  vector[N] m1=X*b;
  for(i in 1:N){
    y1[i]=normal_rng(m1[i],s);
  }
}

normal linear models

n=30
mdl=cmdstan_model('./ex5.stan')

single regression

tb=tibble(x=runif(n,0,9),
          y=rnorm(n,x,1))
f=formula(y~x)
par(mfrow=c(1,1))
plot(tb$x,tb$y)

qplot(data=tb,x,y)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -646.428 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15      -13.4147    0.00105322   0.000473195           1           1       23    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -13.41
##    b[1]      -0.11
##    b[2]       1.00
##    s          0.95
##    y1[1]      5.19
##    y1[2]      9.59
##    y1[3]      4.63
##    y1[4]      5.37
##    y1[5]      4.62
##    y1[6]      6.14
##    y1[7]      8.18
##    y1[8]      6.52
##    y1[9]      3.68
##    y1[10]     1.32
##    y1[11]    -0.74
##    y1[12]     5.57
##    y1[13]     7.96
##    y1[14]     2.44
##    y1[15]     1.89
##    y1[16]     0.68
##    y1[17]     7.27
##    y1[18]     4.85
##    y1[19]     5.32
##    y1[20]     4.03
##    y1[21]     5.97
##    y1[22]     1.60
##    y1[23]     5.44
##    y1[24]     6.08
##    y1[25]     0.90
##    y1[26]     6.35
##    y1[27]     3.06
##    y1[28]     2.60
##    y1[29]     7.14
##    y1[30]     3.89
##    m1[1]      5.92
##    m1[2]      8.31
##    m1[3]      5.58
##    m1[4]      5.05
##    m1[5]      4.43
##    m1[6]      6.15
##    m1[7]      8.33
##    m1[8]      7.61
##    m1[9]      4.53
##    m1[10]     1.00
##    m1[11]     0.07
##    m1[12]     6.96
##    m1[13]     7.58
##    m1[14]     1.84
##    m1[15]     1.62
##    m1[16]     0.48
##    m1[17]     7.64
##    m1[18]     4.03
##    m1[19]     4.00
##    m1[20]     4.95
##    m1[21]     5.15
##    m1[22]     1.18
##    m1[23]     8.66
##    m1[24]     7.98
##    m1[25]     0.89
##    m1[26]     7.73
##    m1[27]     3.24
##    m1[28]     3.40
##    m1[29]     8.27
##    m1[30]     2.88
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -15.00 -14.65 1.29 1.00 -17.71 -13.64 1.01      644      531
##    b[1]    -0.12  -0.11 0.40 0.39  -0.81   0.50 1.01      382      466
##    b[2]     1.00   1.00 0.07 0.07   0.89   1.11 1.01      453      641
##    s        1.02   1.00 0.14 0.13   0.81   1.30 1.01     1115      947
##    y1[1]    5.88   5.87 1.02 1.02   4.18   7.58 1.00     2019     1696
##    y1[2]    8.33   8.31 1.06 1.05   6.59  10.10 1.00     1796     2104
##    y1[3]    5.59   5.57 1.05 1.05   3.88   7.37 1.00     2037     1927
##    y1[4]    5.05   5.06 1.07 0.98   3.28   6.80 1.00     1978     1869
##    y1[5]    4.43   4.44 1.07 1.01   2.67   6.14 1.00     2033     1970
##    y1[6]    6.13   6.12 1.07 1.04   4.38   7.88 1.00     2100     1963
##    y1[7]    8.31   8.32 1.06 1.05   6.57  10.09 1.00     1975     1933
##    y1[8]    7.61   7.62 1.08 1.06   5.85   9.37 1.00     1946     1929
##    y1[9]    4.51   4.52 1.05 1.01   2.83   6.24 1.00     2119     1952
##    y1[10]   1.00   0.99 1.06 1.04  -0.76   2.72 1.00     1534     1968
##    y1[11]   0.08   0.07 1.14 1.10  -1.86   1.92 1.00     1367     1825
##    y1[12]   6.96   6.96 1.03 1.04   5.29   8.68 1.00     1788     1856
##    y1[13]   7.61   7.61 1.06 1.05   5.84   9.31 1.00     2169     1802
##    y1[14]   1.85   1.85 1.09 1.06   0.05   3.69 1.00     1614     1777
##    y1[15]   1.60   1.61 1.11 1.04  -0.21   3.46 1.00     1744     1704
##    y1[16]   0.49   0.48 1.11 1.01  -1.36   2.28 1.00     1433     1965
##    y1[17]   7.65   7.64 1.08 1.06   5.92   9.47 1.00     1838     1772
##    y1[18]   3.99   4.00 1.03 0.98   2.30   5.68 1.00     1860     1774
##    y1[19]   4.05   4.05 1.06 1.02   2.29   5.79 1.00     1895     1933
##    y1[20]   4.96   4.94 1.07 1.06   3.25   6.73 1.00     2030     1950
##    y1[21]   5.18   5.19 1.05 0.99   3.48   6.90 1.00     1945     1849
##    y1[22]   1.17   1.19 1.10 1.04  -0.66   2.91 1.00     1543     1591
##    y1[23]   8.64   8.62 1.05 1.07   6.92  10.33 1.00     1860     1646
##    y1[24]   8.00   7.99 1.06 1.05   6.30   9.80 1.00     1757     1760
##    y1[25]   0.91   0.92 1.08 1.04  -0.84   2.66 1.00     1641     1352
##    y1[26]   7.71   7.73 1.06 1.01   5.94   9.40 1.00     2057     1926
##    y1[27]   3.25   3.28 1.04 1.02   1.50   4.91 1.00     2134     1908
##    y1[28]   3.38   3.44 1.04 1.00   1.64   4.99 1.00     1961     1880
##    y1[29]   8.26   8.27 1.05 1.07   6.51   9.98 1.00     2169     1930
##    y1[30]   2.85   2.87 1.04 1.03   1.15   4.50 1.00     1741     1812
##    m1[1]    5.92   5.91 0.20 0.20   5.61   6.27 1.00     2101     1646
##    m1[2]    8.31   8.31 0.30 0.29   7.83   8.81 1.00     1266     1156
##    m1[3]    5.58   5.57 0.19 0.19   5.27   5.90 1.00     1808     1460
##    m1[4]    5.05   5.05 0.19 0.19   4.75   5.37 1.00     1283     1231
##    m1[5]    4.42   4.42 0.19 0.19   4.12   4.74 1.00      837     1012
##    m1[6]    6.15   6.14 0.20 0.20   5.83   6.50 1.00     2195     1601
##    m1[7]    8.33   8.33 0.30 0.29   7.85   8.83 1.00     1258     1149
##    m1[8]    7.61   7.60 0.26 0.25   7.19   8.05 1.00     1657      994
##    m1[9]    4.53   4.53 0.19 0.18   4.23   4.85 1.00      903     1033
##    m1[10]   0.99   1.01 0.33 0.32   0.42   1.50 1.01      387      478
##    m1[11]   0.07   0.08 0.39 0.38  -0.61   0.67 1.01      383      466
##    m1[12]   6.96   6.96 0.23 0.23   6.59   7.35 1.00     2045     1366
##    m1[13]   7.58   7.58 0.26 0.25   7.17   8.02 1.00     1680      994
##    m1[14]   1.84   1.85 0.28 0.28   1.34   2.27 1.01      400      497
##    m1[15]   1.61   1.63 0.30 0.29   1.10   2.06 1.01      398      496
##    m1[16]   0.47   0.49 0.36 0.36  -0.16   1.03 1.01      384      472
##    m1[17]   7.64   7.64 0.26 0.25   7.22   8.08 1.00     1634      994
##    m1[18]   4.03   4.02 0.20 0.19   3.72   4.35 1.01      669      921
##    m1[19]   4.00   4.00 0.20 0.19   3.69   4.32 1.01      661      895
##    m1[20]   4.95   4.95 0.19 0.19   4.65   5.27 1.00     1200     1200
##    m1[21]   5.15   5.15 0.19 0.19   4.85   5.47 1.00     1371     1232
##    m1[22]   1.18   1.19 0.32 0.32   0.63   1.67 1.01      390      478
##    m1[23]   8.66   8.66 0.32 0.31   8.15   9.19 1.00     1136     1244
##    m1[24]   7.98   7.98 0.28 0.27   7.53   8.45 1.00     1418     1074
##    m1[25]   0.89   0.90 0.34 0.33   0.30   1.40 1.01      386      466
##    m1[26]   7.73   7.73 0.27 0.26   7.30   8.18 1.00     1571     1003
##    m1[27]   3.23   3.24 0.22 0.22   2.87   3.58 1.01      498      657
##    m1[28]   3.39   3.39 0.22 0.21   3.04   3.73 1.01      522      665
##    m1[29]   8.27   8.27 0.30 0.29   7.80   8.77 1.00     1278     1179
##    m1[30]   2.88   2.88 0.23 0.23   2.48   3.25 1.01      448      591

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)

par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

multiple regression

tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
          y=rnorm(n,x1-x2,1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -4204.6 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       25      -14.3992   5.25865e-05   0.000541205           1           1       54    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -14.40
##    b[1]       0.39
##    b[2]       0.93
##    b[3]      -1.05
##    s          0.98
##    y1[1]      3.42
##    y1[2]      6.67
##    y1[3]     -3.01
##    y1[4]      6.78
##    y1[5]     -0.42
##    y1[6]      3.68
##    y1[7]     -4.01
##    y1[8]      1.52
##    y1[9]      4.15
##    y1[10]    -5.76
##    y1[11]    -1.92
##    y1[12]    -2.42
##    y1[13]     1.00
##    y1[14]     5.72
##    y1[15]    -4.48
##    y1[16]     3.54
##    y1[17]     1.75
##    y1[18]    -2.47
##    y1[19]     5.21
##    y1[20]    -2.72
##    y1[21]     0.27
##    y1[22]     1.25
##    y1[23]    -2.80
##    y1[24]     4.40
##    y1[25]    -3.94
##    y1[26]     0.53
##    y1[27]    -0.15
##    y1[28]    -4.72
##    y1[29]     3.05
##    y1[30]    -6.67
##    m1[1]      3.23
##    m1[2]      5.44
##    m1[3]     -2.45
##    m1[4]      7.59
##    m1[5]     -1.23
##    m1[6]      4.17
##    m1[7]     -5.68
##    m1[8]      1.46
##    m1[9]      3.94
##    m1[10]    -5.36
##    m1[11]    -1.49
##    m1[12]    -2.82
##    m1[13]     1.84
##    m1[14]     5.50
##    m1[15]    -3.57
##    m1[16]     3.31
##    m1[17]     1.24
##    m1[18]    -2.68
##    m1[19]     5.61
##    m1[20]    -2.56
##    m1[21]     0.13
##    m1[22]    -0.18
##    m1[23]    -2.46
##    m1[24]     3.92
##    m1[25]    -3.29
##    m1[26]    -0.32
##    m1[27]    -0.87
##    m1[28]    -3.23
##    m1[29]     4.06
##    m1[30]    -6.36
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -16.62 -16.29 1.56 1.38 -19.42 -14.80 1.00      563      791
##    b[1]     0.33   0.34 0.53 0.50  -0.58   1.17 1.01      348      284
##    b[2]     0.94   0.94 0.08 0.08   0.81   1.07 1.01      652      730
##    b[3]    -1.04  -1.04 0.07 0.07  -1.16  -0.93 1.00      799      767
##    s        1.09   1.07 0.16 0.15   0.86   1.38 1.00      983      958
##    y1[1]    3.24   3.27 1.13 1.09   1.29   5.06 1.00     1845     1992
##    y1[2]    5.45   5.44 1.11 1.09   3.61   7.30 1.00     1989     1677
##    y1[3]   -2.53  -2.51 1.20 1.20  -4.45  -0.58 1.00     1532     1712
##    y1[4]    7.56   7.57 1.21 1.12   5.59   9.60 1.00     1884     1854
##    y1[5]   -1.25  -1.22 1.14 1.10  -3.10   0.59 1.00     1575     1842
##    y1[6]    4.20   4.25 1.13 1.11   2.32   6.07 1.00     1621     1689
##    y1[7]   -5.71  -5.67 1.15 1.14  -7.60  -3.91 1.00     1790     1931
##    y1[8]    1.40   1.42 1.17 1.13  -0.52   3.30 1.00     1275     1524
##    y1[9]    3.91   3.88 1.18 1.11   2.08   5.83 1.00     1935     1856
##    y1[10]  -5.35  -5.37 1.16 1.11  -7.25  -3.43 1.00     2018     1812
##    y1[11]  -1.49  -1.50 1.12 1.10  -3.33   0.30 1.00     1901     1697
##    y1[12]  -2.82  -2.82 1.13 1.09  -4.70  -1.02 1.00     1630     1895
##    y1[13]   1.83   1.87 1.18 1.10  -0.22   3.70 1.00     1880     1953
##    y1[14]   5.47   5.44 1.18 1.14   3.63   7.45 1.00     1916     1973
##    y1[15]  -3.58  -3.60 1.21 1.20  -5.54  -1.57 1.00     1753     1974
##    y1[16]   3.31   3.31 1.15 1.12   1.40   5.20 1.00     1551     1572
##    y1[17]   1.27   1.25 1.10 1.07  -0.52   3.08 1.00     2020     2000
##    y1[18]  -2.65  -2.68 1.15 1.13  -4.53  -0.82 1.00     1813     1862
##    y1[19]   5.66   5.70 1.17 1.12   3.70   7.60 1.00     2002     1922
##    y1[20]  -2.52  -2.50 1.14 1.11  -4.35  -0.67 1.00     1988     1880
##    y1[21]   0.09   0.09 1.16 1.10  -1.83   1.94 1.00     1719     1935
##    y1[22]  -0.16  -0.19 1.20 1.20  -2.07   1.79 1.00     1740     1771
##    y1[23]  -2.45  -2.45 1.17 1.19  -4.32  -0.57 1.00     2083     1968
##    y1[24]   3.88   3.89 1.17 1.13   1.98   5.74 1.00     1781     1881
##    y1[25]  -3.34  -3.31 1.13 1.11  -5.20  -1.51 1.00     1829     1948
##    y1[26]  -0.35  -0.33 1.10 1.11  -2.10   1.39 1.00     1839     1404
##    y1[27]  -0.85  -0.87 1.17 1.16  -2.68   1.05 1.00     1771     1894
##    y1[28]  -3.27  -3.27 1.14 1.10  -5.16  -1.42 1.00     2029     1893
##    y1[29]   4.06   4.10 1.12 1.11   2.22   5.90 1.00     1883     1926
##    y1[30]  -6.39  -6.39 1.18 1.15  -8.32  -4.48 1.00     1743     1736
##    m1[1]    3.22   3.22 0.29 0.27   2.74   3.69 1.00      629      739
##    m1[2]    5.44   5.45 0.36 0.35   4.85   6.01 1.00     1600     1019
##    m1[3]   -2.49  -2.49 0.40 0.39  -3.18  -1.87 1.01      386      380
##    m1[4]    7.60   7.60 0.47 0.45   6.81   8.34 1.00     1841     1278
##    m1[5]   -1.27  -1.26 0.34 0.34  -1.84  -0.74 1.01      372      382
##    m1[6]    4.15   4.15 0.35 0.34   3.56   4.70 1.00      566      604
##    m1[7]   -5.69  -5.69 0.38 0.36  -6.32  -5.06 1.00     1609     1392
##    m1[8]    1.42   1.42 0.39 0.36   0.73   2.04 1.01      359      343
##    m1[9]    3.94   3.95 0.30 0.29   3.44   4.43 1.00     1735     1348
##    m1[10]  -5.36  -5.36 0.36 0.35  -5.97  -4.75 1.00     2130     1380
##    m1[11]  -1.49  -1.49 0.24 0.23  -1.91  -1.08 1.00     2153     1440
##    m1[12]  -2.82  -2.82 0.30 0.30  -3.31  -2.32 1.00     1743     1714
##    m1[13]   1.86   1.87 0.34 0.34   1.31   2.39 1.00     1189     1389
##    m1[14]   5.50   5.51 0.37 0.35   4.90   6.08 1.00     1697     1088
##    m1[15]  -3.60  -3.60 0.36 0.36  -4.22  -3.02 1.01      489      587
##    m1[16]   3.28   3.28 0.39 0.38   2.61   3.90 1.01      409      344
##    m1[17]   1.25   1.25 0.23 0.23   0.86   1.62 1.00     1957     1671
##    m1[18]  -2.68  -2.68 0.29 0.28  -3.15  -2.20 1.00     1936     1647
##    m1[19]   5.61   5.61 0.37 0.36   5.00   6.18 1.00     1595     1052
##    m1[20]  -2.54  -2.54 0.38 0.37  -3.16  -1.91 1.00     1047     1245
##    m1[21]   0.09   0.10 0.37 0.36  -0.57   0.68 1.01      354      327
##    m1[22]  -0.15  -0.15 0.46 0.45  -0.90   0.57 1.01      728      915
##    m1[23]  -2.44  -2.45 0.34 0.33  -3.00  -1.88 1.00     1267     1451
##    m1[24]   3.90   3.90 0.32 0.30   3.38   4.41 1.00      661      825
##    m1[25]  -3.31  -3.30 0.29 0.28  -3.79  -2.84 1.00      873     1190
##    m1[26]  -0.34  -0.33 0.24 0.24  -0.73   0.02 1.00      492      555
##    m1[27]  -0.84  -0.84 0.42 0.42  -1.53  -0.17 1.01      801     1081
##    m1[28]  -3.25  -3.24 0.27 0.26  -3.71  -2.80 1.00     1419     1515
##    m1[29]   4.07   4.08 0.32 0.32   3.53   4.59 1.00     1879     1735
##    m1[30]  -6.39  -6.39 0.43 0.42  -7.10  -5.67 1.00      996     1238

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ANOVA

tb=tibble(c=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,(c=='b')*2-(c=='c')*2,1))
f=formula(y~c)
par(mfrow=c(1,1))
boxplot(y~c,tb)

qplot(data=tb,c,y,geom='boxplot')

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -29.1885 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       18      -12.4727   5.90482e-05   6.64483e-05      0.8902      0.8902       21    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -12.47
##    b[1]       0.20
##    b[2]       1.72
##    b[3]      -2.18
##    s          0.92
##    y1[1]     -2.72
##    y1[2]      0.51
##    y1[3]     -0.40
##    y1[4]      0.74
##    y1[5]      0.36
##    y1[6]     -0.08
##    y1[7]      1.91
##    y1[8]      2.80
##    y1[9]      2.74
##    y1[10]     0.24
##    y1[11]     0.06
##    y1[12]     3.02
##    y1[13]    -0.92
##    y1[14]     1.17
##    y1[15]    -2.18
##    y1[16]     1.40
##    y1[17]     0.09
##    y1[18]     2.94
##    y1[19]    -0.20
##    y1[20]     0.74
##    y1[21]     1.68
##    y1[22]    -0.47
##    y1[23]     0.26
##    y1[24]     0.39
##    y1[25]     3.70
##    y1[26]    -1.97
##    y1[27]    -1.66
##    y1[28]     0.91
##    y1[29]     0.47
##    y1[30]     1.82
##    m1[1]     -1.98
##    m1[2]      0.20
##    m1[3]      0.20
##    m1[4]      0.20
##    m1[5]      0.20
##    m1[6]      0.20
##    m1[7]      1.92
##    m1[8]      1.92
##    m1[9]      1.92
##    m1[10]     0.20
##    m1[11]     0.20
##    m1[12]     1.92
##    m1[13]     0.20
##    m1[14]     1.92
##    m1[15]    -1.98
##    m1[16]     1.92
##    m1[17]     0.20
##    m1[18]     1.92
##    m1[19]    -1.98
##    m1[20]     1.92
##    m1[21]     1.92
##    m1[22]     0.20
##    m1[23]     0.20
##    m1[24]    -1.98
##    m1[25]     1.92
##    m1[26]    -1.98
##    m1[27]    -1.98
##    m1[28]     1.92
##    m1[29]     0.20
##    m1[30]     1.92
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -14.71 -14.34 1.53 1.28 -17.67 -12.90 1.00      785     1123
##    b[1]     0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    b[2]     1.73   1.71 0.42 0.40   1.08   2.48 1.00      883      768
##    b[3]    -2.18  -2.18 0.50 0.48  -3.02  -1.38 1.00      682      892
##    s        1.01   1.00 0.15 0.14   0.80   1.27 1.00     2321     1565
##    y1[1]   -1.96  -1.96 1.12 1.10  -3.74  -0.02 1.00     1959     1871
##    y1[2]    0.16   0.15 1.09 1.05  -1.62   1.93 1.00     1816     1751
##    y1[3]    0.20   0.22 1.08 1.06  -1.64   1.99 1.00     1815     2000
##    y1[4]    0.18   0.20 1.06 1.06  -1.54   1.93 1.00     1862     1753
##    y1[5]    0.21   0.24 1.07 1.02  -1.61   1.94 1.00     1777     1638
##    y1[6]    0.20   0.19 1.05 1.00  -1.47   1.95 1.01     1510     1712
##    y1[7]    1.89   1.85 1.12 1.06   0.09   3.73 1.00     1936     1996
##    y1[8]    1.91   1.90 1.06 1.07   0.18   3.64 1.00     1994     1980
##    y1[9]    1.96   1.97 1.06 1.05   0.21   3.69 1.00     1877     1764
##    y1[10]   0.18   0.19 1.04 0.98  -1.60   1.98 1.00     1654     1785
##    y1[11]   0.20   0.20 1.05 1.01  -1.49   1.92 1.00     1598     1955
##    y1[12]   1.90   1.88 1.08 1.03   0.08   3.67 1.00     2085     1875
##    y1[13]   0.20   0.20 1.07 1.08  -1.58   1.90 1.00     2031     1921
##    y1[14]   1.87   1.90 1.10 1.03   0.08   3.60 1.00     1994     1805
##    y1[15]  -1.98  -1.99 1.08 1.05  -3.82  -0.22 1.00     1754     1615
##    y1[16]   1.92   1.90 1.08 1.05   0.17   3.68 1.00     2115     1554
##    y1[17]   0.18   0.19 1.09 1.04  -1.61   2.00 1.00     1484     1807
##    y1[18]   1.95   1.97 1.09 1.02   0.10   3.78 1.00     1755     1929
##    y1[19]  -2.02  -2.02 1.14 1.10  -3.92  -0.15 1.00     1651     1958
##    y1[20]   1.92   1.91 1.08 1.03   0.16   3.71 1.00     1610     1795
##    y1[21]   1.91   1.92 1.08 1.02   0.13   3.64 1.00     2075     1957
##    y1[22]   0.20   0.22 1.07 1.06  -1.56   1.94 1.00     1790     1804
##    y1[23]   0.20   0.20 1.07 1.06  -1.56   1.98 1.00     1672     1643
##    y1[24]  -2.00  -2.01 1.12 1.11  -3.83  -0.18 1.00     1915     1772
##    y1[25]   1.94   1.94 1.05 1.02   0.17   3.60 1.00     1681     1908
##    y1[26]  -1.97  -1.94 1.11 1.12  -3.83  -0.14 1.00     1677     1881
##    y1[27]  -1.99  -1.98 1.07 1.07  -3.75  -0.30 1.00     1818     1690
##    y1[28]   1.93   1.93 1.04 1.01   0.19   3.68 1.00     1970     1953
##    y1[29]   0.21   0.21 1.06 1.07  -1.54   1.94 1.00     1516     1456
##    y1[30]   1.96   1.93 1.08 1.06   0.18   3.75 1.00     1753     1871
##    m1[1]   -1.98  -1.98 0.41 0.39  -2.65  -1.32 1.00     1092     1331
##    m1[2]    0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[3]    0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[4]    0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[5]    0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[6]    0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[7]    1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[8]    1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[9]    1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[10]   0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[11]   0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[12]   1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[13]   0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[14]   1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[15]  -1.98  -1.98 0.41 0.39  -2.65  -1.32 1.00     1092     1331
##    m1[16]   1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[17]   0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[18]   1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[19]  -1.98  -1.98 0.41 0.39  -2.65  -1.32 1.00     1092     1331
##    m1[20]   1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[21]   1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[22]   0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[23]   0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[24]  -1.98  -1.98 0.41 0.39  -2.65  -1.32 1.00     1092     1331
##    m1[25]   1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[26]  -1.98  -1.98 0.41 0.39  -2.65  -1.32 1.00     1092     1331
##    m1[27]  -1.98  -1.98 0.41 0.39  -2.65  -1.32 1.00     1092     1331
##    m1[28]   1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272
##    m1[29]   0.19   0.20 0.30 0.31  -0.32   0.67 1.00      867     1175
##    m1[30]   1.92   1.92 0.30 0.28   1.46   2.41 1.00     1508     1272

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')

qplot(data=tb,y,y1,shape=c,size=I(2),xlab='obs.',ylab='prd.')

plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')

qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(data=tb,1:n,y,col=c)+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ANCOVA

tb=tibble(x=runif(n,0,9),c=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,x+(c=='b')*2-(c=='c')*2,1))

f=formula(y~x+c)
par(mfrow=c(1,1))
#plot(tb$x1,tb$y,pch=tb$c1)

lv=c(0,1,2)
plot(tb$x,tb$y,pch=lv[factor(tb$c)])

qplot(data=tb,x,y,shape=c,size=I(2))

plot(tb$x,tb$y,col=1+lv[factor(tb$c)])

qplot(data=tb,x,y,col=c)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -152.843 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       40      -8.25775   5.84488e-05    0.00125024           1           1       49    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      -8.26
##    b[1]      -0.08
##    b[2]       0.95
##    b[3]       2.86
##    b[4]      -1.84
##    s          0.80
##    y1[1]      4.76
##    y1[2]     -1.93
##    y1[3]      0.72
##    y1[4]      0.93
##    y1[5]      5.97
##    y1[6]      8.09
##    y1[7]     10.46
##    y1[8]      7.13
##    y1[9]      5.33
##    y1[10]     6.06
##    y1[11]    -0.12
##    y1[12]    -1.08
##    y1[13]     6.56
##    y1[14]     5.25
##    y1[15]     2.12
##    y1[16]     6.44
##    y1[17]     4.81
##    y1[18]     0.49
##    y1[19]     5.53
##    y1[20]     4.02
##    y1[21]     8.40
##    y1[22]     0.32
##    y1[23]     4.53
##    y1[24]     4.66
##    y1[25]     4.79
##    y1[26]     4.16
##    y1[27]     7.74
##    y1[28]     8.27
##    y1[29]     5.39
##    y1[30]     2.88
##    m1[1]      4.30
##    m1[2]     -1.68
##    m1[3]      1.11
##    m1[4]      0.57
##    m1[5]      4.59
##    m1[6]      8.26
##    m1[7]     11.17
##    m1[8]      8.04
##    m1[9]      5.40
##    m1[10]     6.54
##    m1[11]    -0.12
##    m1[12]    -0.64
##    m1[13]     6.19
##    m1[14]     5.64
##    m1[15]     1.70
##    m1[16]     6.49
##    m1[17]     4.00
##    m1[18]     1.01
##    m1[19]     4.75
##    m1[20]     2.92
##    m1[21]     7.76
##    m1[22]     1.05
##    m1[23]     4.55
##    m1[24]     4.51
##    m1[25]     4.23
##    m1[26]     3.76
##    m1[27]     7.58
##    m1[28]     8.59
##    m1[29]     5.39
##    m1[30]     3.28
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -11.19 -10.84 1.72 1.55 -14.42 -9.09 1.01      751     1048
##    b[1]    -0.09  -0.09 0.34 0.33  -0.64  0.49 1.00      925      949
##    b[2]     0.95   0.95 0.06 0.06   0.85  1.04 1.00     1606     1661
##    b[3]     2.86   2.85 0.37 0.35   2.28  3.47 1.01      814     1121
##    b[4]    -1.81  -1.81 0.45 0.44  -2.56 -1.07 1.01      799      857
##    s        0.91   0.89 0.13 0.12   0.71  1.14 1.00     2143     1743
##    y1[1]    4.30   4.30 0.94 0.93   2.72  5.79 1.00     1959     1990
##    y1[2]   -1.66  -1.68 1.05 1.02  -3.35  0.09 1.00     1692     1814
##    y1[3]    1.09   1.11 0.98 0.92  -0.56  2.73 1.00     1511     1868
##    y1[4]    0.53   0.54 0.96 0.91  -1.04  2.06 1.00     1938     1985
##    y1[5]    4.61   4.62 0.96 0.92   3.00  6.22 1.00     1837     1903
##    y1[6]    8.24   8.23 0.96 0.93   6.74  9.83 1.00     1770     2031
##    y1[7]   11.18  11.19 1.00 1.00   9.58 12.80 1.00     1951     1917
##    y1[8]    8.01   8.03 0.96 0.92   6.43  9.54 1.00     1850     1894
##    y1[9]    5.35   5.35 0.98 0.97   3.73  6.92 1.00     1865     1902
##    y1[10]   6.53   6.52 1.01 1.00   4.90  8.14 1.00     1933     1958
##    y1[11]  -0.11  -0.10 1.01 0.96  -1.75  1.50 1.00     1501     1866
##    y1[12]  -0.64  -0.64 1.00 0.97  -2.25  0.98 1.00     1797     1855
##    y1[13]   6.19   6.17 0.94 0.90   4.69  7.77 1.00     1642     1728
##    y1[14]   5.67   5.66 1.02 0.98   3.95  7.32 1.00     2042     1954
##    y1[15]   1.70   1.70 0.95 0.91   0.14  3.23 1.00     1759     1893
##    y1[16]   6.46   6.46 0.97 0.95   4.86  8.01 1.00     2145     1963
##    y1[17]   4.00   3.99 0.94 0.89   2.47  5.52 1.00     1891     1842
##    y1[18]   0.98   0.95 1.00 1.02  -0.63  2.63 1.00     1733     1584
##    y1[19]   4.77   4.80 0.95 0.89   3.24  6.30 1.00     1939     1807
##    y1[20]   2.93   2.90 0.98 0.93   1.37  4.59 1.00     1874     1855
##    y1[21]   7.72   7.72 0.95 0.93   6.13  9.32 1.00     1923     1969
##    y1[22]   1.04   1.04 0.93 0.92  -0.51  2.56 1.00     1846     1906
##    y1[23]   4.52   4.55 0.98 0.91   2.92  6.14 1.00     2052     1919
##    y1[24]   4.51   4.53 0.97 0.95   2.94  6.10 1.00     1733     1863
##    y1[25]   4.19   4.20 0.94 0.94   2.66  5.71 1.00     1929     1930
##    y1[26]   3.73   3.71 0.93 0.93   2.18  5.26 1.00     1985     1933
##    y1[27]   7.56   7.57 0.98 0.94   5.97  9.21 1.00     1809     1848
##    y1[28]   8.60   8.61 0.99 0.97   6.99 10.18 1.00     1960     1920
##    y1[29]   5.36   5.37 0.94 0.94   3.78  6.86 1.00     1909     1940
##    y1[30]   3.32   3.34 0.99 0.97   1.63  4.94 1.00     1803     1766
##    m1[1]    4.29   4.28 0.30 0.29   3.80  4.80 1.00     1471     1122
##    m1[2]   -1.67  -1.67 0.46 0.43  -2.42 -0.89 1.01     1077     1000
##    m1[3]    1.09   1.09 0.30 0.29   0.63  1.59 1.00      875     1144
##    m1[4]    0.55   0.55 0.32 0.31   0.05  1.08 1.00      891     1031
##    m1[5]    4.58   4.58 0.25 0.25   4.16  4.98 1.01      992     1322
##    m1[6]    8.25   8.25 0.31 0.29   7.74  8.75 1.00     1228     1410
##    m1[7]   11.15  11.15 0.42 0.39  10.46 11.84 1.00     1303     1267
##    m1[8]    8.03   8.02 0.30 0.28   7.54  8.53 1.00     1227     1411
##    m1[9]    5.38   5.38 0.28 0.27   4.94  5.85 1.00     1372     1288
##    m1[10]   6.55   6.56 0.46 0.45   5.78  7.27 1.01     1405     1407
##    m1[11]  -0.11  -0.11 0.42 0.38  -0.80  0.58 1.01     1056      885
##    m1[12]  -0.63  -0.63 0.43 0.40  -1.34  0.09 1.01     1060      826
##    m1[13]   6.17   6.17 0.28 0.26   5.73  6.63 1.00     1265     1306
##    m1[14]   5.66   5.67 0.43 0.42   4.93  6.32 1.01     1332     1292
##    m1[15]   1.69   1.69 0.28 0.26   1.25  2.15 1.00      852     1034
##    m1[16]   6.48   6.49 0.30 0.30   5.96  6.96 1.01     1287     1254
##    m1[17]   3.99   4.00 0.25 0.24   3.58  4.38 1.01      926     1203
##    m1[18]   0.99   0.99 0.30 0.29   0.52  1.50 1.00      881     1114
##    m1[19]   4.74   4.75 0.26 0.25   4.32  5.15 1.01     1017     1271
##    m1[20]   2.90   2.89 0.34 0.33   2.33  3.49 1.00     1570     1211
##    m1[21]   7.75   7.76 0.36 0.35   7.15  8.32 1.01     1503     1321
##    m1[22]   1.03   1.03 0.30 0.29   0.56  1.53 1.00      878     1117
##    m1[23]   4.53   4.52 0.29 0.28   4.06  5.03 1.00     1452     1192
##    m1[24]   4.50   4.49 0.30 0.28   4.02  5.00 1.00     1454     1192
##    m1[25]   4.22   4.22 0.25 0.24   3.80  4.61 1.01      949     1251
##    m1[26]   3.75   3.75 0.25 0.24   3.34  4.14 1.01      904     1225
##    m1[27]   7.57   7.59 0.35 0.34   6.99  8.13 1.01     1464     1326
##    m1[28]   8.58   8.58 0.31 0.30   8.05  9.10 1.00     1238     1331
##    m1[29]   5.38   5.39 0.27 0.26   4.92  5.81 1.01     1109     1282
##    m1[30]   3.29   3.31 0.39 0.37   2.65  3.91 1.01     1150      994

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)

plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')+
  geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(data=tb,1:n,y,col=c)+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ex6

generalized linear regression

poisson regression

objective variable [0,Infinity), integer. variance of error is near to mean  
(normal linear regression, correlation is near to 1,-1, variance of error is constant)  

# of samples is N,  
log li=sum(bj*xji),j=[0,K],i=[1,N]  
yi~Po(li)  
 or  
li=sum(bj xji),j=[0,k]  
yi~Po(exp li)  

when xj -> xj +1, y -> y* exp bj   

ex6-1.stan

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] l=X*b;
  y~poisson_log(l);
}
generated quantities{
  array[N] int y1;
  vector[N] l1=X*b;
  for(i in 1:N){
    y1[i]=poisson_log_rng(l1[i]);
  }
}
n=30
tb=tibble(x=runif(n,-1,1),c=sample(c('a','b'),n,replace=T),
          y=rpois(n,exp(x+(c=='b')*0.5)))
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

glm(f,tb,family='poisson')
## 
## Call:  glm(formula = f, family = "poisson", data = tb)
## 
## Coefficients:
## (Intercept)            x           cb  
##      0.0658       1.0269       0.7378  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       48.5 
## Residual Deviance: 26.7  AIC: 90.2
X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -61.439 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        8      -13.0485   0.000132489   0.000119888           1           1       11    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -13.05
##    b[1]       0.07
##    b[2]       1.03
##    b[3]       0.74
##    y1[1]      0.00
##    y1[2]      0.00
##    y1[3]      2.00
##    y1[4]      5.00
##    y1[5]      0.00
##    y1[6]      2.00
##    y1[7]      0.00
##    y1[8]      7.00
##    y1[9]      2.00
##    y1[10]     2.00
##    y1[11]     1.00
##    y1[12]     1.00
##    y1[13]     0.00
##    y1[14]     1.00
##    y1[15]     1.00
##    y1[16]     0.00
##    y1[17]     0.00
##    y1[18]     1.00
##    y1[19]     2.00
##    y1[20]     4.00
##    y1[21]     3.00
##    y1[22]     1.00
##    y1[23]     1.00
##    y1[24]     0.00
##    y1[25]     6.00
##    y1[26]     1.00
##    y1[27]     3.00
##    y1[28]     2.00
##    y1[29]     7.00
##    y1[30]     1.00
##    l1[1]     -0.09
##    l1[2]     -0.17
##    l1[3]      0.27
##    l1[4]      1.66
##    l1[5]     -0.23
##    l1[6]      0.98
##    l1[7]      0.62
##    l1[8]      1.04
##    l1[9]      0.89
##    l1[10]     0.10
##    l1[11]     0.16
##    l1[12]     0.03
##    l1[13]    -0.64
##    l1[14]     0.07
##    l1[15]     0.23
##    l1[16]    -0.14
##    l1[17]    -0.18
##    l1[18]     0.53
##    l1[19]     0.33
##    l1[20]     1.63
##    l1[21]     0.04
##    l1[22]    -0.90
##    l1[23]     0.23
##    l1[24]    -0.77
##    l1[25]     0.85
##    l1[26]     0.55
##    l1[27]     0.60
##    l1[28]     0.64
##    l1[29]     1.41
##    l1[30]    -0.06
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='l1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -14.53 -14.30 1.11 0.98 -16.54 -13.25 1.00     1024     1083
##    b[1]     0.04   0.04 0.21 0.22  -0.30   0.38 1.00     1110     1283
##    b[2]     1.03   1.02 0.25 0.26   0.63   1.44 1.00     1289     1447
##    b[3]     0.74   0.74 0.28 0.29   0.27   1.18 1.00     1073     1312
##    y1[1]    0.96   1.00 1.05 1.48   0.00   3.00 1.00     1755     1590
##    y1[2]    0.82   1.00 0.93 1.48   0.00   3.00 1.00     2041     1929
##    y1[3]    1.34   1.00 1.15 1.48   0.00   3.00 1.00     1937     1960
##    y1[4]    5.20   5.00 2.63 2.97   1.00  10.00 1.00     1732     1962
##    y1[5]    0.81   1.00 0.92 1.48   0.00   3.00 1.00     1888     1862
##    y1[6]    2.67   2.00 1.71 1.48   0.00   6.00 1.00     1933     2055
##    y1[7]    1.78   2.00 1.40 1.48   0.00   4.00 1.00     2060     2048
##    y1[8]    2.84   3.00 1.76 1.48   0.00   6.00 1.00     1923     1968
##    y1[9]    2.33   2.00 1.61 1.48   0.00   5.00 1.00     1891     1881
##    y1[10]   1.11   1.00 1.09 1.48   0.00   3.00 1.00     2036     1987
##    y1[11]   1.20   1.00 1.15 1.48   0.00   3.00 1.00     1925     1893
##    y1[12]   1.01   1.00 1.04 1.48   0.00   3.00 1.00     1812     1577
##    y1[13]   0.57   0.00 0.80 0.00   0.00   2.00 1.00     1842     2025
##    y1[14]   1.12   1.00 1.08 1.48   0.00   3.00 1.00     1863     1843
##    y1[15]   1.23   1.00 1.12 1.48   0.00   3.00 1.00     1884     1909
##    y1[16]   0.86   1.00 0.95 1.48   0.00   3.00 1.00     2002     1982
##    y1[17]   0.90   1.00 0.99 1.48   0.00   3.00 1.00     2004     1813
##    y1[18]   1.69   1.00 1.34 1.48   0.00   4.00 1.00     2030     1915
##    y1[19]   1.36   1.00 1.24 1.48   0.00   4.00 1.00     2005     1990
##    y1[20]   5.10   5.00 2.56 2.97   1.00  10.00 1.00     1902     2084
##    y1[21]   1.00   1.00 1.00 1.48   0.00   3.00 1.00     2086     2021
##    y1[22]   0.40   0.00 0.66 0.00   0.00   2.00 1.00     1791     1817
##    y1[23]   1.23   1.00 1.11 1.48   0.00   3.00 1.00     1890     1916
##    y1[24]   0.45   0.00 0.71 0.00   0.00   2.00 1.00     1905     1948
##    y1[25]   2.29   2.00 1.61 1.48   0.00   5.00 1.00     2111     1993
##    y1[26]   1.71   2.00 1.38 1.48   0.00   4.00 1.00     2126     1936
##    y1[27]   1.88   2.00 1.46 1.48   0.00   5.00 1.00     1892     1949
##    y1[28]   1.86   2.00 1.41 1.48   0.00   4.00 1.00     1887     1984
##    y1[29]   4.19   4.00 2.30 2.97   1.00   8.00 1.00     1616     1904
##    y1[30]   0.90   1.00 0.97 1.48   0.00   3.00 1.00     2042     1935
##    l1[1]   -0.11  -0.10 0.34 0.33  -0.68   0.42 1.00     1467     1537
##    l1[2]   -0.20  -0.20 0.24 0.25  -0.58   0.20 1.00     1078     1096
##    l1[3]    0.24   0.24 0.19 0.20  -0.07   0.57 1.00     1194     1230
##    l1[4]    1.63   1.63 0.25 0.26   1.20   2.04 1.00     1452     1564
##    l1[5]   -0.26  -0.26 0.25 0.26  -0.65   0.15 1.00     1076     1150
##    l1[6]    0.95   0.95 0.23 0.24   0.56   1.33 1.00     1536     1334
##    l1[7]    0.59   0.59 0.19 0.20   0.26   0.90 1.00     1397     1377
##    l1[8]    1.01   1.01 0.24 0.25   0.61   1.40 1.00     1543     1359
##    l1[9]    0.86   0.86 0.22 0.23   0.49   1.22 1.00     1517     1475
##    l1[10]   0.07   0.08 0.30 0.30  -0.44   0.56 1.00     1509     1533
##    l1[11]   0.13   0.14 0.29 0.29  -0.37   0.60 1.00     1525     1585
##    l1[12]   0.00   0.00 0.22 0.22  -0.34   0.35 1.00     1101     1248
##    l1[13]  -0.67  -0.67 0.33 0.33  -1.19  -0.12 1.00     1082     1367
##    l1[14]   0.04   0.05 0.31 0.30  -0.48   0.53 1.00     1499     1533
##    l1[15]   0.20   0.20 0.20 0.20  -0.13   0.52 1.00     1166     1264
##    l1[16]  -0.17  -0.17 0.24 0.24  -0.54   0.22 1.00     1079     1122
##    l1[17]  -0.21  -0.20 0.36 0.35  -0.81   0.35 1.00     1452     1479
##    l1[18]   0.50   0.51 0.24 0.24   0.10   0.88 1.00     1662     1639
##    l1[19]   0.30   0.30 0.19 0.20  -0.01   0.62 1.00     1225     1230
##    l1[20]   1.60   1.61 0.25 0.25   1.18   2.01 1.00     1464     1566
##    l1[21]   0.01   0.01 0.21 0.22  -0.34   0.36 1.00     1102     1248
##    l1[22]  -0.93  -0.93 0.38 0.39  -1.54  -0.29 1.00     1099     1378
##    l1[23]   0.20   0.20 0.20 0.20  -0.12   0.53 1.00     1166     1280
##    l1[24]  -0.80  -0.80 0.35 0.36  -1.36  -0.21 1.00     1091     1363
##    l1[25]   0.82   0.82 0.21 0.22   0.46   1.17 1.00     1506     1545
##    l1[26]   0.52   0.52 0.19 0.20   0.19   0.83 1.00     1358     1377
##    l1[27]   0.57   0.58 0.23 0.23   0.19   0.94 1.00     1698     1654
##    l1[28]   0.61   0.62 0.22 0.23   0.23   0.97 1.00     1716     1602
##    l1[29]   1.38   1.39 0.22 0.23   1.01   1.75 1.00     1566     1383
##    l1[30]  -0.09  -0.09 0.23 0.23  -0.44   0.28 1.00     1086     1217

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##     0 1 2 3 4 5
##   0 2 6 0 0 0 0
##   1 1 5 1 0 0 0
##   2 0 5 1 1 1 0
##   3 0 0 4 0 0 0
##   4 0 0 0 0 0 1
##   6 0 0 1 0 0 1
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     0 1 2 3 4 5
##   0 2 4 0 0 0 0
##   1 1 4 1 0 0 0
##   2 0 2 1 1 0 0
##   3 0 0 2 0 0 0
##   4 0 0 0 0 0 0
##   6 0 0 1 0 0 0
## 
## , ,  = b
## 
##    
##     0 1 2 3 4 5
##   0 0 2 0 0 0 0
##   1 0 1 0 0 0 0
##   2 0 3 0 0 1 0
##   3 0 0 2 0 0 0
##   4 0 0 0 0 0 1
##   6 0 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##     0 1 2 3 4
##   0 5 3 0 0 0
##   1 4 3 0 0 0
##   2 1 4 2 1 0
##   3 0 3 1 0 0
##   4 0 0 0 0 1
##   6 0 0 1 0 1
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     0 1 2 3 4
##   0 4 2 0 0 0
##   1 4 2 0 0 0
##   2 0 2 2 0 0
##   3 0 1 1 0 0
##   4 0 0 0 0 0
##   6 0 0 1 0 0
## 
## , ,  = b
## 
##    map
##     0 1 2 3 4
##   0 1 1 0 0 0
##   1 0 1 0 0 0
##   2 1 2 0 1 0
##   3 0 2 0 0 0
##   4 0 0 0 0 1
##   6 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

logistic regression   

# of samples is N,  
objective variable 0/1 binary  
  
probability of incident pi[0,1]  
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)  

yi~Ber(pi), 0/1 binary  

odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident  
odds ratio(x0,x1)=odds(x1)/odd(x0)  
  
when xj -> xj +1, odds ratio -> odds ratio *exp bj  

ex6-2.stan

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] z=X*b;
  y~bernoulli_logit(z);
}
generated quantities{
  array[N] int y1;
  vector[N] z1=X*b;
  for(i in 1:N){
    y1[i]=bernoulli_rng(inv_logit(z1[i]));
  }
}
n=30
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,1,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)

f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

glm(f,tb,family='binomial') # it can caluculte when all trials are once
## 
## Call:  glm(formula = f, family = "binomial", data = tb)
## 
## Coefficients:
## (Intercept)            x           cb  
##       0.018        4.034        2.779  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       32.6 
## Residual Deviance: 20.7  AIC: 26.7
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -28.9854 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11      -10.3708   0.000328813   2.39543e-05           1           1       15    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -10.37
##    b[1]       0.02
##    b[2]       4.03
##    b[3]       2.78
##    y1[1]      0.00
##    y1[2]      1.00
##    y1[3]      1.00
##    y1[4]      0.00
##    y1[5]      1.00
##    y1[6]      1.00
##    y1[7]      0.00
##    y1[8]      1.00
##    y1[9]      1.00
##    y1[10]     0.00
##    y1[11]     1.00
##    y1[12]     1.00
##    y1[13]     0.00
##    y1[14]     1.00
##    y1[15]     1.00
##    y1[16]     0.00
##    y1[17]     1.00
##    y1[18]     1.00
##    y1[19]     1.00
##    y1[20]     1.00
##    y1[21]     1.00
##    y1[22]     1.00
##    y1[23]     1.00
##    y1[24]     1.00
##    y1[25]     1.00
##    y1[26]     0.00
##    y1[27]     1.00
##    y1[28]     0.00
##    y1[29]     1.00
##    y1[30]     1.00
##    z1[1]      1.68
##    z1[2]      1.09
##    z1[3]      2.46
##    z1[4]     -2.51
##    z1[5]      0.61
##    z1[6]      4.20
##    z1[7]      0.61
##    z1[8]      3.71
##    z1[9]      2.52
##    z1[10]     1.36
##    z1[11]    -0.69
##    z1[12]     0.18
##    z1[13]    -1.02
##    z1[14]     1.75
##    z1[15]     3.27
##    z1[16]     1.72
##    z1[17]     6.08
##    z1[18]     5.77
##    z1[19]     1.09
##    z1[20]     2.85
##    z1[21]     4.88
##    z1[22]     3.86
##    z1[23]     6.13
##    z1[24]     6.31
##    z1[25]     4.12
##    z1[26]     0.26
##    z1[27]     0.80
##    z1[28]    -2.05
##    z1[29]     2.03
##    z1[30]     1.01
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='z1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -12.03 -11.67 1.35 1.07 -14.84 -10.58 1.00      666      807
##    b[1]    -0.13  -0.08 0.96 0.94  -1.81   1.35 1.01      710      759
##    b[2]     5.28   5.00 2.19 2.07   2.19   9.30 1.00      414      544
##    b[3]     3.72   3.50 1.95 1.83   0.95   7.34 1.01      397      544
##    y1[1]    0.86   1.00 0.35 0.00   0.00   1.00 1.00     1908       NA
##    y1[2]    0.75   1.00 0.43 0.00   0.00   1.00 1.00     1934       NA
##    y1[3]    0.92   1.00 0.27 0.00   0.00   1.00 1.00     1980       NA
##    y1[4]    0.09   0.00 0.28 0.00   0.00   1.00 1.00     1131       NA
##    y1[5]    0.65   1.00 0.48 0.00   0.00   1.00 1.00     1829       NA
##    y1[6]    0.98   1.00 0.15 0.00   1.00   1.00 1.00     1849       NA
##    y1[7]    0.66   1.00 0.47 0.00   0.00   1.00 1.00     1761       NA
##    y1[8]    0.98   1.00 0.15 0.00   1.00   1.00 1.00     1772       NA
##    y1[9]    0.93   1.00 0.26 0.00   0.00   1.00 1.00     1723       NA
##    y1[10]   0.81   1.00 0.39 0.00   0.00   1.00 1.00     1881       NA
##    y1[11]   0.31   0.00 0.46 0.00   0.00   1.00 1.00     1924       NA
##    y1[12]   0.55   1.00 0.50 0.00   0.00   1.00 1.00     2037       NA
##    y1[13]   0.24   0.00 0.42 0.00   0.00   1.00 1.00     1332       NA
##    y1[14]   0.86   1.00 0.35 0.00   0.00   1.00 1.00     1783       NA
##    y1[15]   0.96   1.00 0.20 0.00   1.00   1.00 1.00     1938       NA
##    y1[16]   0.86   1.00 0.34 0.00   0.00   1.00 1.00     1931       NA
##    y1[17]   1.00   1.00 0.07 0.00   1.00   1.00 1.00     1689       NA
##    y1[18]   1.00   1.00 0.06 0.00   1.00   1.00 1.00     2030       NA
##    y1[19]   0.77   1.00 0.42 0.00   0.00   1.00 1.00     1961       NA
##    y1[20]   0.95   1.00 0.22 0.00   1.00   1.00 1.00     1835       NA
##    y1[21]   0.99   1.00 0.10 0.00   1.00   1.00 1.00     1717       NA
##    y1[22]   0.98   1.00 0.15 0.00   1.00   1.00 1.00     1712       NA
##    y1[23]   0.99   1.00 0.07 0.00   1.00   1.00 1.00     1192       NA
##    y1[24]   1.00   1.00 0.07 0.00   1.00   1.00 1.00     1736       NA
##    y1[25]   0.98   1.00 0.13 0.00   1.00   1.00 1.00     1932       NA
##    y1[26]   0.56   1.00 0.50 0.00   0.00   1.00 1.00     2165       NA
##    y1[27]   0.70   1.00 0.46 0.00   0.00   1.00 1.00     2141       NA
##    y1[28]   0.13   0.00 0.33 0.00   0.00   1.00 1.00     1548       NA
##    y1[29]   0.88   1.00 0.32 0.00   0.00   1.00 1.00     2003       NA
##    y1[30]   0.76   1.00 0.43 0.00   0.00   1.00 1.00     1941       NA
##    z1[1]    2.04   1.95 1.01 0.95   0.57   3.78 1.00     1098     1137
##    z1[2]    1.27   1.21 0.89 0.87  -0.07   2.78 1.00     1201     1203
##    z1[3]    3.07   2.97 1.27 1.23   1.24   5.29 1.00      771      872
##    z1[4]   -3.45  -3.20 1.98 1.87  -7.05  -0.65 1.01      411      544
##    z1[5]    0.64   0.62 0.88 0.86  -0.75   2.09 1.00     1018     1117
##    z1[6]    5.43   5.15 2.07 1.90   2.58   9.28 1.00      392      559
##    z1[7]    0.73   0.70 0.84 0.83  -0.65   2.15 1.00     1878     1404
##    z1[8]    4.70   4.53 1.83 1.78   2.10   7.97 1.00      578      761
##    z1[9]    3.22   3.07 1.29 1.20   1.41   5.55 1.00      455      651
##    z1[10]   1.70   1.64 0.91 0.89   0.29   3.29 1.00      806      979
##    z1[11]  -0.98  -0.94 1.14 1.08  -2.93   0.77 1.00     1055     1055
##    z1[12]   0.16   0.14 0.89 0.86  -1.28   1.61 1.00     1891     1443
##    z1[13]  -1.49  -1.33 1.31 1.20  -3.83   0.43 1.01      482      607
##    z1[14]   2.14   2.05 1.03 0.97   0.64   3.91 1.00     1064     1137
##    z1[15]   4.21   3.99 1.62 1.49   1.96   7.20 1.00      410      569
##    z1[16]   2.18   2.08 1.00 0.97   0.70   3.95 1.00      611      764
##    z1[17]   7.89   7.48 3.03 2.78   3.76  13.46 1.00      381      541
##    z1[18]   7.48   7.12 2.87 2.64   3.55  12.74 1.00      382      541
##    z1[19]   1.27   1.21 0.89 0.87  -0.07   2.78 1.00     1201     1203
##    z1[20]   3.57   3.44 1.43 1.39   1.54   6.05 1.00      683      711
##    z1[21]   6.32   6.00 2.41 2.20   3.00  10.84 1.00      385      556
##    z1[22]   4.89   4.72 1.90 1.85   2.20   8.30 1.00      567      706
##    z1[23]   7.95   7.52 3.05 2.80   3.77  13.56 1.00      381      541
##    z1[24]   8.19   7.74 3.15 2.87   3.88  13.99 1.00      381      541
##    z1[25]   5.32   5.05 2.03 1.86   2.52   9.08 1.00      393      559
##    z1[26]   0.27   0.25 0.88 0.84  -1.14   1.72 1.00     1953     1519
##    z1[27]   0.97   0.94 0.84 0.82  -0.41   2.39 1.00     1561     1354
##    z1[28]  -2.84  -2.63 1.76 1.66  -5.96  -0.34 1.01      423      533
##    z1[29]   2.51   2.42 1.12 1.07   0.90   4.40 1.00      920      984
##    z1[30]   1.25   1.21 0.86 0.83  -0.11   2.71 1.00     1185     1251

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##      0  1
##   0  3  4
##   1  1 22
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##      0  1
##   0  3  1
##   1  0  9
## 
## , ,  = b
## 
##    
##      0  1
##   0  0  3
##   1  1 13
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##      0  1
##   0  3  4
##   1  1 22
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##      0  1
##   0  3  1
##   1  0  9
## 
## , ,  = b
## 
##    map
##      0  1
##   0  0  3
##   1  1 13
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

multi logistic regression  

# of samples is N,  
# of trials of a sample i is mi,  
objective variable [0,n], integer  
  
probability of incident pi[0,1]  
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)  

yi~B(mi, pi), # of acutual incident  

odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident  
odds ratio(x0,x1)=odds(x1)/odd(x0)  
  
when xj -> xj +1, odds ratio -> odds ratio *exp bj  

ex6-3.stan

data{
  int N;
  int K;
  array[N] int m;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] z=X*b;
  y~binomial_logit(m,z);
}
generated quantities{
  array[N] int y1;
  vector[N] z1=X*b;
  for(i in 1:N){
    y1[i]=binomial_rng(m[i],inv_logit(z1[i]));
  }
}
n=30
m=floor(runif(n,1,10)) # trials are varying (1,10)
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,m,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)

f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X,m=m)

mdl=cmdstan_model('./ex6-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -180.062 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        9      -114.455    0.00117686    0.00149793           1           1       11    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__    -114.46
##    b[1]       0.12
##    b[2]       0.43
##    b[3]       0.38
##    y1[1]      4.00
##    y1[2]      3.00
##    y1[3]      2.00
##    y1[4]      4.00
##    y1[5]      0.00
##    y1[6]      2.00
##    y1[7]      6.00
##    y1[8]      5.00
##    y1[9]      4.00
##    y1[10]     4.00
##    y1[11]     1.00
##    y1[12]     1.00
##    y1[13]     4.00
##    y1[14]     2.00
##    y1[15]     4.00
##    y1[16]     4.00
##    y1[17]     4.00
##    y1[18]     6.00
##    y1[19]     4.00
##    y1[20]     7.00
##    y1[21]     4.00
##    y1[22]     5.00
##    y1[23]     2.00
##    y1[24]     2.00
##    y1[25]     8.00
##    y1[26]     3.00
##    y1[27]     2.00
##    y1[28]     4.00
##    y1[29]     0.00
##    y1[30]     1.00
##    z1[1]      0.43
##    z1[2]      0.09
##    z1[3]      0.15
##    z1[4]      0.36
##    z1[5]     -0.06
##    z1[6]     -0.27
##    z1[7]     -0.23
##    z1[8]      0.77
##    z1[9]     -0.11
##    z1[10]     0.56
##    z1[11]     0.04
##    z1[12]     0.59
##    z1[13]     0.87
##    z1[14]     0.25
##    z1[15]     0.08
##    z1[16]     0.34
##    z1[17]     0.91
##    z1[18]     0.23
##    z1[19]     0.45
##    z1[20]     0.33
##    z1[21]     0.42
##    z1[22]     0.40
##    z1[23]     0.47
##    z1[24]     0.71
##    z1[25]     0.87
##    z1[26]     0.53
##    z1[27]     0.42
##    z1[28]     0.61
##    z1[29]    -0.04
##    z1[30]     0.22
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='z1',ch=F)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -116.01 -115.68 1.33 1.00 -118.54 -114.62 1.00      826     1232
##    b[1]      0.12    0.13 0.22 0.21   -0.24    0.48 1.00     1264     1200
##    b[2]      0.44    0.45 0.28 0.29   -0.02    0.89 1.00     1353     1256
##    b[3]      0.38    0.38 0.32 0.32   -0.12    0.91 1.00      972     1007
##    y1[1]     3.61    4.00 1.24 1.48    2.00    6.00 1.00     2026       NA
##    y1[2]     3.06    3.00 1.27 1.48    1.00    5.00 1.00     1929     1998
##    y1[3]     1.06    1.00 0.71 1.48    0.00    2.00 1.00     1956       NA
##    y1[4]     2.34    2.00 0.99 1.48    1.00    4.00 1.00     1853       NA
##    y1[5]     2.94    3.00 1.26 1.48    1.00    5.00 1.00     1849     1973
##    y1[6]     2.56    3.00 1.30 1.48    0.95    5.00 1.00     2094     2053
##    y1[7]     3.49    3.00 1.52 1.48    1.00    6.00 1.00     1806     1836
##    y1[8]     5.45    6.00 1.41 1.48    3.00    8.00 1.00     1836       NA
##    y1[9]     2.83    3.00 1.27 1.48    1.00    5.00 1.00     1857     1890
##    y1[10]    3.18    3.00 1.12 1.48    1.00    5.00 1.00     1913       NA
##    y1[11]    1.56    2.00 0.88 1.48    0.00    3.00 1.00     1985       NA
##    y1[12]    1.28    1.00 0.68 1.48    0.00    2.00 1.00     1853       NA
##    y1[13]    4.22    4.00 1.17 1.48    2.00    6.00 1.00     1961       NA
##    y1[14]    4.53    5.00 1.44 1.48    2.00    7.00 1.00     1785     1893
##    y1[15]    2.08    2.00 1.05 1.48    0.00    4.00 1.00     1863       NA
##    y1[16]    4.73    5.00 1.46 1.48    2.00    7.00 1.00     1926     1953
##    y1[17]    3.58    4.00 1.06 1.48    2.00    5.00 1.00     1774       NA
##    y1[18]    5.04    5.00 1.66 1.48    2.00    8.00 1.00     1724     1743
##    y1[19]    5.52    6.00 1.54 1.48    3.00    8.00 1.00     1990     1796
##    y1[20]    5.22    5.00 1.57 1.48    3.00    8.00 1.00     1905     1683
##    y1[21]    3.04    3.00 1.11 1.48    1.00    5.00 1.00     1758       NA
##    y1[22]    4.19    4.00 1.38 1.48    2.00    6.00 1.00     2076     1956
##    y1[23]    2.49    3.00 0.98 1.48    1.00    4.00 1.00     1924       NA
##    y1[24]    2.69    3.00 0.96 1.48    1.00    4.00 1.00     1713       NA
##    y1[25]    6.35    6.00 1.42 1.48    4.00    9.00 1.00     1817       NA
##    y1[26]    3.78    4.00 1.25 1.48    2.00    6.00 1.00     1833       NA
##    y1[27]    1.21    1.00 0.69 1.48    0.00    2.00 1.00     1608       NA
##    y1[28]    5.90    6.00 1.48 1.48    3.00    8.00 1.00     1784     1784
##    y1[29]    1.48    1.00 0.87 1.48    0.00    3.00 1.00     2088       NA
##    y1[30]    1.66    2.00 0.89 1.48    0.00    3.00 1.00     1947       NA
##    z1[1]     0.44    0.43 0.27 0.27    0.00    0.89 1.00     1369     1538
##    z1[2]     0.09    0.10 0.22 0.21   -0.27    0.45 1.00     1259     1333
##    z1[3]     0.16    0.15 0.36 0.36   -0.42    0.74 1.00     1427     1259
##    z1[4]     0.37    0.37 0.27 0.28   -0.07    0.81 1.00     1527     1618
##    z1[5]    -0.06   -0.05 0.26 0.25   -0.50    0.37 1.00     1252     1127
##    z1[6]    -0.27   -0.27 0.35 0.34   -0.85    0.31 1.00     1279     1233
##    z1[7]    -0.23   -0.23 0.33 0.32   -0.77    0.32 1.00     1271     1165
##    z1[8]     0.79    0.78 0.27 0.26    0.37    1.24 1.00     1555     1744
##    z1[9]    -0.11   -0.11 0.28 0.27   -0.58    0.34 1.00     1255     1138
##    z1[10]    0.57    0.57 0.23 0.23    0.20    0.96 1.00     1578     1616
##    z1[11]    0.04    0.05 0.23 0.22   -0.34    0.42 1.00     1252     1129
##    z1[12]    0.60    0.60 0.23 0.23    0.23    0.98 1.00     1572     1692
##    z1[13]    0.89    0.88 0.30 0.31    0.40    1.40 1.00     1496     1672
##    z1[14]    0.26    0.26 0.22 0.21   -0.10    0.62 1.00     1295     1264
##    z1[15]    0.08    0.08 0.40 0.40   -0.56    0.73 1.00     1410     1276
##    z1[16]    0.35    0.35 0.24 0.24   -0.04    0.75 1.00     1333     1426
##    z1[17]    0.93    0.92 0.32 0.33    0.42    1.48 1.00     1492     1586
##    z1[18]    0.24    0.24 0.32 0.33   -0.29    0.76 1.00     1453     1273
##    z1[19]    0.46    0.45 0.28 0.28    0.00    0.93 1.00     1375     1602
##    z1[20]    0.34    0.33 0.24 0.23   -0.05    0.73 1.00     1329     1403
##    z1[21]    0.43    0.42 0.27 0.26   -0.01    0.87 1.00     1366     1591
##    z1[22]    0.41    0.41 0.26 0.26   -0.01    0.83 1.00     1556     1526
##    z1[23]    0.48    0.47 0.24 0.24    0.08    0.87 1.00     1604     1457
##    z1[24]    0.73    0.72 0.25 0.25    0.33    1.15 1.00     1570     1781
##    z1[25]    0.89    0.88 0.31 0.31    0.40    1.41 1.00     1494     1672
##    z1[26]    0.54    0.53 0.32 0.31    0.04    1.07 1.00     1392     1598
##    z1[27]    0.43    0.43 0.25 0.25    0.02    0.84 1.00     1569     1520
##    z1[28]    0.63    0.62 0.23 0.23    0.25    1.02 1.00     1567     1739
##    z1[29]   -0.04   -0.03 0.25 0.24   -0.46    0.38 1.00     1253     1062
##    z1[30]    0.22    0.22 0.22 0.21   -0.13    0.57 1.00     1283     1197

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##     1 2 3 4 5 6
##   1 2 2 0 0 0 0
##   2 2 1 2 1 0 0
##   3 0 0 5 1 1 0
##   4 0 1 2 2 0 0
##   5 0 0 0 1 2 2
##   6 0 0 0 0 1 1
##   7 0 0 0 0 0 1
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     1 2 3 4 5 6
##   1 0 2 0 0 0 0
##   2 1 0 1 0 0 0
##   3 0 0 3 0 1 0
##   4 0 0 2 2 0 0
##   5 0 0 0 0 2 0
##   6 0 0 0 0 0 1
##   7 0 0 0 0 0 0
## 
## , ,  = b
## 
##    
##     1 2 3 4 5 6
##   1 2 0 0 0 0 0
##   2 1 1 1 1 0 0
##   3 0 0 2 1 0 0
##   4 0 1 0 0 0 0
##   5 0 0 0 1 0 2
##   6 0 0 0 0 1 0
##   7 0 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##     1 2 3 4 5 6 7
##   1 2 2 0 0 0 0 0
##   2 2 1 2 1 0 0 0
##   3 0 0 5 2 0 0 0
##   4 0 1 2 2 0 0 0
##   5 0 0 0 0 3 2 0
##   6 0 0 0 0 1 1 0
##   7 0 0 0 0 0 0 1
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     1 2 3 4 5 6 7
##   1 0 2 0 0 0 0 0
##   2 1 0 1 0 0 0 0
##   3 0 0 3 1 0 0 0
##   4 0 0 2 2 0 0 0
##   5 0 0 0 0 2 0 0
##   6 0 0 0 0 0 1 0
##   7 0 0 0 0 0 0 0
## 
## , ,  = b
## 
##    map
##     1 2 3 4 5 6 7
##   1 2 0 0 0 0 0 0
##   2 1 1 1 1 0 0 0
##   3 0 0 2 1 0 0 0
##   4 0 1 0 0 0 0 0
##   5 0 0 0 0 1 2 0
##   6 0 0 0 0 1 0 0
##   7 0 0 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

ex7 interaction of variables

n=50
mdl=cmdstan_model('./ex5.stan')

tb=tibble(x=runif(n,-3,3),
          ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
          y=rnorm(n,x+(ca=='b')+(cb=='b')-(ca=='b')*(cb=='b'),1))

grid.arrange(qplot(data=tb,x,y,col=ca),
             qplot(data=tb,x,y,col=cb),ncol=2)

fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
mle

mcmc=goMCMC(mdl,data)

seeMCMC(mcmc,exc='m',ch=F)

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

grid.arrange(
  qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
grid.arrange(
  qplot(data=tb,1:n,y,col=ca)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  qplot(data=tb,1:n,y,col=cb)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  ncol=2
)
}


f0=formula(y~x+ca)
f1=formula(y~x+ca+cb)
f2=formula(y~x+ca*cb)

fn(f0)
## Initial log joint probability = -613.308 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       14      -30.2492    0.00026474   0.000991645      0.9285      0.9285       15    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -32.27 -31.94 1.50 1.24 -35.31 -30.48 1.00      918     1472
##    b[1]     0.34   0.34 0.23 0.22  -0.06   0.70 1.01      755      948
##    b[2]     0.90   0.90 0.11 0.11   0.72   1.08 1.00     2517     1595
##    b[3]     0.68   0.68 0.34 0.34   0.12   1.26 1.01      608      707
##    s        1.18   1.17 0.13 0.12   0.99   1.41 1.00     2424     1683
##    y1[1]   -0.76  -0.73 1.22 1.22  -2.84   1.18 1.00     1904     1933
##    y1[2]    0.15   0.11 1.24 1.27  -1.81   2.20 1.00     1743     1599
##    y1[3]    0.43   0.42 1.21 1.19  -1.57   2.46 1.00     1789     2011
##    y1[4]    0.58   0.58 1.18 1.18  -1.35   2.47 1.00     1803     2054
##    y1[5]   -0.29  -0.30 1.22 1.21  -2.24   1.72 1.00     1771     1704
##    y1[6]   -0.16  -0.16 1.22 1.17  -2.22   1.85 1.00     2031     2016
##    y1[7]    0.18   0.20 1.26 1.22  -1.89   2.22 1.00     1783     1856
##    y1[8]    3.38   3.39 1.24 1.25   1.34   5.53 1.00     2173     2055
##    y1[9]   -1.88  -1.90 1.24 1.18  -3.89   0.14 1.00     2059     1970
##    y1[10]   0.35   0.36 1.21 1.20  -1.65   2.33 1.00     1914     1893
##    y1[11]   0.47   0.47 1.22 1.20  -1.50   2.49 1.00     1869     1792
##    y1[12]   1.16   1.17 1.18 1.13  -0.82   3.15 1.00     2006     1753
##    y1[13]   2.08   2.09 1.22 1.22   0.08   4.07 1.00     2063     2129
##    y1[14]   1.77   1.79 1.23 1.19  -0.25   3.74 1.00     1986     1966
##    y1[15]   1.91   1.89 1.24 1.17  -0.12   3.97 1.00     1854     1865
##    y1[16]   3.06   3.03 1.23 1.20   1.09   5.08 1.00     1953     1928
##    y1[17]   3.00   2.98 1.24 1.25   1.03   5.04 1.00     1991     1852
##    y1[18]   1.58   1.58 1.24 1.20  -0.43   3.66 1.00     2020     1657
##    y1[19]  -0.33  -0.33 1.24 1.21  -2.36   1.72 1.00     1944     2057
##    y1[20]   3.53   3.51 1.24 1.20   1.39   5.57 1.00     2117     1955
##    y1[21]  -1.59  -1.64 1.19 1.18  -3.52   0.37 1.00     1820     1973
##    y1[22]   0.23   0.25 1.17 1.13  -1.67   2.15 1.00     1941     1712
##    y1[23]   1.95   1.96 1.21 1.19  -0.06   3.90 1.00     1577     1527
##    y1[24]   2.61   2.62 1.24 1.27   0.56   4.66 1.00     1558     1694
##    y1[25]   1.09   1.10 1.24 1.17  -1.00   3.08 1.00     2028     1914
##    y1[26]  -1.38  -1.39 1.20 1.25  -3.29   0.58 1.00     2065     1967
##    y1[27]  -0.10  -0.03 1.23 1.21  -2.16   1.88 1.00     1994     1831
##    y1[28]   2.13   2.12 1.19 1.18   0.19   4.04 1.00     2023     2013
##    y1[29]   2.31   2.32 1.21 1.21   0.39   4.38 1.00     1912     1779
##    y1[30]   1.70   1.70 1.19 1.20  -0.33   3.63 1.00     1960     1971
##    y1[31]  -1.11  -1.11 1.23 1.16  -3.16   0.93 1.00     1857     1811
##    y1[32]  -1.50  -1.51 1.23 1.20  -3.53   0.53 1.00     2004     1933
##    y1[33]   2.55   2.52 1.28 1.24   0.44   4.67 1.00     1746     1841
##    y1[34]   0.21   0.23 1.21 1.24  -1.82   2.25 1.00     1915     1630
##    y1[35]   2.67   2.66 1.19 1.15   0.66   4.63 1.00     1771     1822
##    y1[36]  -1.35  -1.35 1.23 1.18  -3.44   0.65 1.00     1847     1781
##    y1[37]  -0.28  -0.31 1.23 1.26  -2.32   1.66 1.00     1827     1672
##    y1[38]   0.35   0.35 1.19 1.16  -1.60   2.30 1.00     1828     1855
##    y1[39]  -2.18  -2.19 1.25 1.18  -4.21  -0.16 1.00     1782     1623
##    y1[40]  -0.29  -0.27 1.19 1.21  -2.24   1.62 1.00     1863     1858
##    y1[41]   0.89   0.94 1.22 1.18  -1.09   2.91 1.00     1968     1892
##    y1[42]  -0.95  -0.94 1.26 1.24  -3.06   1.11 1.00     1862     1933
##    y1[43]   3.10   3.10 1.23 1.21   1.10   5.13 1.00     2015     1894
##    y1[44]   1.42   1.44 1.21 1.17  -0.60   3.44 1.00     1905     1928
##    y1[45]   1.91   1.93 1.18 1.17  -0.02   3.90 1.00     2188     1664
##    y1[46]  -0.73  -0.76 1.24 1.21  -2.75   1.32 1.00     2122     1839
##    y1[47]   0.00  -0.03 1.22 1.16  -1.97   2.12 1.00     1727     1868
##    y1[48]   2.04   2.01 1.21 1.15   0.08   4.06 1.00     1876     1954
##    y1[49]   1.81   1.82 1.19 1.14  -0.22   3.67 1.00     1832     1884
##    y1[50]   2.92   2.92 1.20 1.22   0.94   4.85 1.00     2095     1973
##    m1[1]   -0.71  -0.71 0.26 0.26  -1.14  -0.28 1.00      996     1185
##    m1[2]    0.11   0.12 0.23 0.22  -0.28   0.48 1.01      778      950
##    m1[3]    0.41   0.42 0.23 0.22   0.02   0.78 1.01      751      988
##    m1[4]    0.56   0.57 0.23 0.23   0.16   0.93 1.01      750      901
##    m1[5]   -0.34  -0.34 0.33 0.33  -0.86   0.19 1.00     1189     1411
##    m1[6]   -0.14  -0.13 0.23 0.23  -0.53   0.24 1.01      824      977
##    m1[7]    0.23   0.23 0.29 0.29  -0.22   0.70 1.00     1081     1309
##    m1[8]    3.33   3.33 0.34 0.33   2.77   3.89 1.00     1762     1532
##    m1[9]   -1.88  -1.88 0.35 0.36  -2.45  -1.33 1.00     1525     1540
##    m1[10]   0.35   0.35 0.28 0.28  -0.09   0.81 1.00     1063     1327
##    m1[11]   0.53   0.53 0.23 0.23   0.13   0.89 1.01      750      913
##    m1[12]   1.17   1.17 0.25 0.25   0.76   1.57 1.00     1022     1448
##    m1[13]   2.09   2.09 0.26 0.25   1.68   2.52 1.00     1221     1593
##    m1[14]   1.76   1.76 0.25 0.24   1.35   2.17 1.00     1115     1426
##    m1[15]   1.88   1.88 0.30 0.30   1.39   2.37 1.00      936     1162
##    m1[16]   3.06   3.05 0.32 0.31   2.53   3.59 1.00     1633     1522
##    m1[17]   2.97   2.96 0.31 0.30   2.46   3.49 1.00     1587     1558
##    m1[18]   1.60   1.60 0.25 0.24   1.19   2.01 1.00     1078     1432
##    m1[19]  -0.38  -0.38 0.33 0.33  -0.90   0.17 1.00     1196     1390
##    m1[20]   3.56   3.55 0.36 0.35   2.97   4.15 1.00     1860     1529
##    m1[21]  -1.53  -1.53 0.32 0.32  -2.04  -1.02 1.00     1359     1466
##    m1[22]   0.23   0.24 0.23 0.22  -0.16   0.60 1.01      763      968
##    m1[23]   1.93   1.94 0.31 0.31   1.44   2.43 1.00      949     1191
##    m1[24]   2.58   2.58 0.37 0.37   1.99   3.18 1.00     1104     1328
##    m1[25]   1.13   1.14 0.25 0.25   0.72   1.55 1.00      796      972
##    m1[26]  -1.39  -1.39 0.31 0.31  -1.89  -0.90 1.00     1292     1451
##    m1[27]  -0.07  -0.07 0.23 0.23  -0.46   0.30 1.01      811      987
##    m1[28]   2.13   2.13 0.26 0.25   1.71   2.56 1.00     1234     1619
##    m1[29]   2.32   2.31 0.27 0.26   1.88   2.76 1.00     1308     1549
##    m1[30]   1.72   1.73 0.29 0.29   1.25   2.19 1.00      901     1162
##    m1[31]  -1.13  -1.13 0.29 0.29  -1.61  -0.66 1.00     1172     1350
##    m1[32]  -1.51  -1.51 0.44 0.44  -2.21  -0.80 1.00     1440     1496
##    m1[33]   2.48   2.48 0.36 0.36   1.90   3.06 1.00     1082     1227
##    m1[34]   0.23   0.23 0.29 0.29  -0.22   0.70 1.00     1081     1309
##    m1[35]   2.67   2.67 0.29 0.28   2.20   3.15 1.00     1461     1545
##    m1[36]  -1.34  -1.34 0.42 0.42  -2.01  -0.66 1.00     1407     1454
##    m1[37]  -0.27  -0.26 0.24 0.24  -0.66   0.12 1.01      855     1030
##    m1[38]   0.35   0.35 0.23 0.22  -0.05   0.71 1.01      754      948
##    m1[39]  -2.16  -2.17 0.38 0.39  -2.78  -1.56 1.00     1651     1495
##    m1[40]  -0.30  -0.29 0.24 0.24  -0.70   0.09 1.01      864     1022
##    m1[41]   0.94   0.94 0.24 0.24   0.53   1.33 1.00      771      924
##    m1[42]  -0.91  -0.92 0.38 0.38  -1.52  -0.30 1.00     1312     1442
##    m1[43]   3.13   3.13 0.32 0.32   2.59   3.67 1.00     1664     1484
##    m1[44]   1.44   1.44 0.25 0.24   1.03   1.85 1.00     1050     1422
##    m1[45]   1.90   1.90 0.31 0.30   1.40   2.39 1.00      941     1149
##    m1[46]  -0.76  -0.75 0.26 0.26  -1.19  -0.32 1.00     1014     1215
##    m1[47]   0.04   0.05 0.23 0.22  -0.35   0.41 1.01      789      977
##    m1[48]   2.03   2.03 0.26 0.25   1.62   2.45 1.00     1198     1593
##    m1[49]   1.82   1.82 0.25 0.24   1.41   2.23 1.00     1131     1427
##    m1[50]   2.93   2.93 0.31 0.30   2.43   3.45 1.00     1570     1558

fn(f1)
## Initial log joint probability = -85.83 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16      -26.1273   0.000562122    0.00136671           1           1       27    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -28.81 -28.44 1.69 1.45 -32.08 -26.73 1.00      878      927
##    b[1]    -0.04  -0.03 0.26 0.25  -0.48   0.40 1.00     1038     1177
##    b[2]     0.92   0.92 0.10 0.10   0.75   1.09 1.00     2389     1402
##    b[3]     0.46   0.45 0.32 0.31  -0.06   1.00 1.00      844     1257
##    b[4]     0.91   0.91 0.32 0.31   0.38   1.42 1.00     1099     1149
##    s        1.10   1.08 0.12 0.12   0.92   1.31 1.00     2119     1463
##    y1[1]   -0.16  -0.16 1.14 1.11  -2.04   1.68 1.00     1757     1967
##    y1[2]   -0.27  -0.32 1.16 1.08  -2.12   1.66 1.00     1881     1475
##    y1[3]    0.98   1.02 1.14 1.13  -0.87   2.84 1.00     2000     1810
##    y1[4]    0.16   0.17 1.12 1.09  -1.69   1.99 1.00     1844     1895
##    y1[5]   -0.05  -0.04 1.18 1.16  -1.97   1.88 1.00     2007     1846
##    y1[6]   -0.55  -0.58 1.14 1.14  -2.38   1.30 1.00     1942     1932
##    y1[7]    0.52   0.51 1.09 1.07  -1.20   2.34 1.00     2033     1972
##    y1[8]    2.81   2.82 1.16 1.16   0.88   4.74 1.00     1914     1921
##    y1[9]   -1.38  -1.38 1.15 1.15  -3.29   0.46 1.00     2010     1925
##    y1[10]  -0.29  -0.30 1.18 1.17  -2.27   1.66 1.00     1972     1678
##    y1[11]   0.15   0.16 1.12 1.09  -1.64   1.96 1.00     1886     1805
##    y1[12]   1.50   1.52 1.11 1.16  -0.38   3.28 1.00     2099     2013
##    y1[13]   2.46   2.48 1.15 1.11   0.57   4.35 1.00     1939     1715
##    y1[14]   1.18   1.16 1.12 1.05  -0.68   2.99 1.00     1962     1973
##    y1[15]   2.48   2.46 1.18 1.15   0.58   4.44 1.00     1831     1819
##    y1[16]   3.43   3.43 1.15 1.14   1.59   5.38 1.00     1888     1874
##    y1[17]   2.44   2.46 1.13 1.10   0.56   4.24 1.00     2041     1933
##    y1[18]   0.97   1.02 1.16 1.14  -0.93   2.86 1.01     1876     1886
##    y1[19]  -0.13  -0.11 1.17 1.14  -2.05   1.76 1.00     1862     1722
##    y1[20]   3.91   3.89 1.15 1.08   2.08   5.85 1.00     2010     1798
##    y1[21]  -1.92  -1.90 1.15 1.11  -3.86  -0.01 1.00     1803     1695
##    y1[22]   0.84   0.84 1.12 1.10  -1.08   2.72 1.00     1870     1783
##    y1[23]   1.59   1.55 1.15 1.12  -0.26   3.53 1.00     1860     1890
##    y1[24]   3.14   3.14 1.19 1.20   1.17   5.01 1.00     1891     1863
##    y1[25]   1.68   1.69 1.15 1.13  -0.22   3.57 1.00     1907     1833
##    y1[26]  -0.89  -0.89 1.14 1.15  -2.73   0.90 1.00     1899     1737
##    y1[27]  -0.45  -0.45 1.15 1.13  -2.30   1.45 1.00     1853     1882
##    y1[28]   2.48   2.48 1.15 1.14   0.59   4.31 1.00     1926     1918
##    y1[29]   2.65   2.63 1.14 1.10   0.81   4.57 1.00     1908     1768
##    y1[30]   2.25   2.26 1.15 1.13   0.41   4.17 1.00     2060     1695
##    y1[31]  -1.59  -1.57 1.14 1.13  -3.51   0.29 1.00     1918     1961
##    y1[32]  -2.19  -2.21 1.18 1.15  -4.16  -0.28 1.00     1974     1784
##    y1[33]   2.17   2.20 1.17 1.14   0.21   4.04 1.00     1922     1932
##    y1[34]   0.52   0.52 1.12 1.14  -1.30   2.34 1.00     1818     1834
##    y1[35]   2.12   2.09 1.16 1.16   0.29   4.01 1.00     1845     1818
##    y1[36]  -1.06  -1.07 1.19 1.16  -2.99   0.88 1.00     1873     1882
##    y1[37]  -0.70  -0.69 1.14 1.14  -2.55   1.16 1.00     1939     1850
##    y1[38]  -0.02   0.00 1.14 1.12  -1.94   1.76 1.00     1804     1956
##    y1[39]  -2.61  -2.61 1.18 1.13  -4.55  -0.69 1.00     1942     1813
##    y1[40]   0.21   0.22 1.13 1.10  -1.69   2.01 1.00     1770     1921
##    y1[41]   0.62   0.60 1.11 1.07  -1.22   2.43 1.00     1858     1869
##    y1[42]  -0.67  -0.69 1.17 1.16  -2.59   1.23 1.00     1884     1937
##    y1[43]   2.59   2.57 1.14 1.15   0.75   4.52 1.00     1925     2013
##    y1[44]   1.76   1.76 1.13 1.10  -0.04   3.66 1.00     1920     2014
##    y1[45]   2.47   2.43 1.15 1.09   0.60   4.37 1.00     1838     1854
##    y1[46]  -1.15  -1.15 1.16 1.16  -3.05   0.74 1.00     2093     2013
##    y1[47]  -0.33  -0.33 1.16 1.19  -2.23   1.51 1.00     1935     1888
##    y1[48]   2.38   2.39 1.15 1.13   0.51   4.28 1.00     2090     1956
##    y1[49]   2.16   2.17 1.12 1.09   0.33   3.98 1.00     2003     2030
##    y1[50]   2.33   2.31 1.17 1.12   0.41   4.30 1.00     1715     1714
##    m1[1]   -0.20  -0.20 0.31 0.31  -0.69   0.31 1.00     1286     1428
##    m1[2]   -0.26  -0.26 0.26 0.25  -0.72   0.17 1.00     1048     1195
##    m1[3]    0.95   0.96 0.29 0.28   0.48   1.44 1.00     1163     1321
##    m1[4]    0.19   0.19 0.26 0.25  -0.24   0.63 1.00     1037     1162
##    m1[5]   -0.06  -0.06 0.32 0.33  -0.59   0.47 1.00     1652     1528
##    m1[6]   -0.52  -0.52 0.27 0.25  -0.98  -0.08 1.00     1077     1340
##    m1[7]    0.52   0.53 0.29 0.30   0.05   0.99 1.00     1566     1557
##    m1[8]    2.79   2.79 0.36 0.36   2.18   3.38 1.00     1577     1474
##    m1[9]   -1.39  -1.40 0.37 0.37  -2.01  -0.79 1.00     1608     1491
##    m1[10]  -0.26  -0.26 0.33 0.33  -0.81   0.26 1.00     1217     1454
##    m1[11]   0.16   0.16 0.26 0.25  -0.27   0.60 1.00     1037     1162
##    m1[12]   1.49   1.49 0.26 0.26   1.06   1.91 1.00     1543     1513
##    m1[13]   2.43   2.43 0.27 0.28   2.00   2.87 1.00     1780     1521
##    m1[14]   1.18   1.18 0.30 0.30   0.69   1.65 1.00     1217     1367
##    m1[15]   2.45   2.44 0.35 0.33   1.90   3.03 1.00     1327     1382
##    m1[16]   3.42   3.42 0.32 0.33   2.91   3.93 1.00     2143     1484
##    m1[17]   2.41   2.42 0.34 0.34   1.85   2.97 1.00     1476     1451
##    m1[18]   1.02   1.02 0.30 0.30   0.53   1.50 1.00     1192     1399
##    m1[19]  -0.10  -0.09 0.33 0.33  -0.64   0.44 1.00     1659     1650
##    m1[20]   3.93   3.93 0.36 0.36   3.35   4.51 1.00     2335     1424
##    m1[21]  -1.94  -1.94 0.34 0.33  -2.49  -1.37 1.00     1487     1550
##    m1[22]   0.77   0.77 0.29 0.28   0.30   1.26 1.00     1167     1235
##    m1[23]   1.60   1.60 0.32 0.32   1.07   2.12 1.00     1198     1247
##    m1[24]   3.17   3.16 0.40 0.39   2.52   3.83 1.00     1463     1526
##    m1[25]   1.69   1.69 0.31 0.30   1.20   2.21 1.00     1206     1247
##    m1[26]  -0.89  -0.90 0.34 0.34  -1.45  -0.33 1.00     1463     1510
##    m1[27]  -0.46  -0.46 0.27 0.25  -0.91  -0.02 1.00     1068     1251
##    m1[28]   2.47   2.47 0.27 0.28   2.03   2.91 1.00     1790     1454
##    m1[29]   2.66   2.66 0.28 0.28   2.21   3.12 1.00     1859     1454
##    m1[30]   2.29   2.28 0.34 0.32   1.75   2.86 1.00     1296     1370
##    m1[31]  -1.54  -1.54 0.31 0.30  -2.05  -1.01 1.00     1326     1501
##    m1[32]  -2.16  -2.15 0.46 0.46  -2.92  -1.43 1.00     1589     1548
##    m1[33]   2.15   2.16 0.36 0.35   1.58   2.76 1.00     1310     1346
##    m1[34]   0.53   0.53 0.29 0.30   0.05   1.00 1.00     1565     1557
##    m1[35]   2.11   2.11 0.33 0.33   1.58   2.63 1.00     1401     1482
##    m1[36]  -1.08  -1.07 0.40 0.41  -1.74  -0.42 1.00     1812     1608
##    m1[37]  -0.65  -0.65 0.27 0.25  -1.11  -0.20 1.00     1098     1360
##    m1[38]  -0.02  -0.02 0.26 0.25  -0.47   0.41 1.00     1037     1162
##    m1[39]  -2.59  -2.60 0.39 0.39  -3.23  -1.96 1.00     1746     1578
##    m1[40]   0.22   0.22 0.29 0.29  -0.26   0.71 1.00     1209     1340
##    m1[41]   0.58   0.57 0.27 0.26   0.12   1.02 1.00     1059     1224
##    m1[42]  -0.64  -0.64 0.37 0.37  -1.25  -0.04 1.00     1756     1660
##    m1[43]   2.58   2.58 0.35 0.35   1.99   3.14 1.00     1520     1473
##    m1[44]   1.76   1.77 0.26 0.27   1.34   2.19 1.00     1583     1545
##    m1[45]   2.47   2.46 0.35 0.34   1.92   3.06 1.00     1331     1382
##    m1[46]  -1.15  -1.15 0.29 0.28  -1.64  -0.66 1.00     1209     1457
##    m1[47]  -0.33  -0.34 0.26 0.25  -0.79   0.10 1.00     1054     1226
##    m1[48]   2.37   2.37 0.27 0.27   1.94   2.80 1.00     1755     1541
##    m1[49]   2.15   2.16 0.26 0.27   1.73   2.58 1.00     1687     1552
##    m1[50]   2.38   2.39 0.34 0.34   1.81   2.93 1.00     1467     1388

fn(f2)
## Initial log joint probability = -82.744 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       25      -24.0366   0.000100239    0.00122831           1           1       30    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -27.29 -26.94 1.91 1.63 -31.04 -24.93 1.01      638      850
##    b[1]    -0.30  -0.29 0.29 0.28  -0.81   0.18 1.00      600      862
##    b[2]     0.88   0.87 0.10 0.10   0.72   1.03 1.00     1814     1588
##    b[3]     1.19   1.19 0.49 0.48   0.40   1.98 1.01      473      824
##    b[4]     1.49   1.47 0.44 0.43   0.78   2.20 1.01      461      670
##    b[5]    -1.27  -1.23 0.68 0.68  -2.40  -0.19 1.01      388      457
##    s        1.07   1.06 0.12 0.11   0.90   1.28 1.00     1720     1524
##    y1[1]    0.17   0.17 1.11 1.10  -1.67   1.99 1.00     2006     1909
##    y1[2]   -0.49  -0.48 1.12 1.09  -2.33   1.31 1.00     1601     2013
##    y1[3]    1.23   1.24 1.12 1.08  -0.56   3.02 1.00     1692     1997
##    y1[4]   -0.10  -0.10 1.11 1.06  -1.94   1.75 1.00     1870     1973
##    y1[5]   -0.19  -0.19 1.12 1.05  -2.06   1.70 1.00     2056     2059
##    y1[6]   -0.77  -0.77 1.11 1.06  -2.60   1.01 1.00     1735     1954
##    y1[7]    0.29   0.30 1.13 1.15  -1.56   2.21 1.00     2024     1783
##    y1[8]    3.15   3.15 1.13 1.09   1.27   4.93 1.00     1889     1820
##    y1[9]   -0.97  -0.97 1.18 1.15  -2.91   0.96 1.00     1669     1407
##    y1[10]   0.23   0.23 1.16 1.18  -1.68   2.14 1.00     1740     1954
##    y1[11]  -0.09  -0.11 1.09 1.06  -1.84   1.76 1.00     1735     2001
##    y1[12]   1.25   1.27 1.08 1.08  -0.48   3.03 1.00     2032     1876
##    y1[13]   2.16   2.14 1.11 1.12   0.35   4.01 1.00     1774     1884
##    y1[14]   1.62   1.62 1.14 1.11  -0.15   3.50 1.00     1903     1761
##    y1[15]   2.68   2.69 1.16 1.13   0.74   4.62 1.00     1970     1971
##    y1[16]   3.10   3.12 1.13 1.13   1.24   4.91 1.00     1862     1739
##    y1[17]   2.78   2.78 1.11 1.06   0.93   4.63 1.00     1823     2055
##    y1[18]   1.45   1.43 1.15 1.10  -0.38   3.32 1.00     1679     1916
##    y1[19]  -0.28  -0.27 1.13 1.09  -2.15   1.52 1.00     2051     2127
##    y1[20]   3.56   3.56 1.12 1.09   1.69   5.38 1.00     1977     1954
##    y1[21]  -2.14  -2.16 1.14 1.10  -3.98  -0.26 1.00     1916     1840
##    y1[22]   1.08   1.08 1.11 1.09  -0.78   2.90 1.00     1925     1968
##    y1[23]   1.26   1.25 1.13 1.13  -0.62   3.10 1.00     1820     1787
##    y1[24]   3.39   3.40 1.15 1.12   1.40   5.23 1.00     1962     1876
##    y1[25]   1.94   1.96 1.10 1.07   0.12   3.75 1.00     1661     1703
##    y1[26]  -0.50  -0.50 1.13 1.12  -2.38   1.35 1.00     1711     1786
##    y1[27]  -0.71  -0.72 1.12 1.11  -2.50   1.16 1.00     1690     1627
##    y1[28]   2.19   2.23 1.10 1.06   0.34   3.99 1.00     1977     1895
##    y1[29]   2.34   2.37 1.10 1.15   0.54   4.10 1.00     1785     1885
##    y1[30]   2.53   2.55 1.10 1.06   0.70   4.32 1.00     1867     1700
##    y1[31]  -1.73  -1.72 1.11 1.09  -3.59   0.07 1.00     2087     2030
##    y1[32]  -1.57  -1.60 1.18 1.17  -3.51   0.39 1.00     1847     2009
##    y1[33]   1.78   1.82 1.14 1.14  -0.14   3.60 1.00     1836     1745
##    y1[34]   0.33   0.35 1.13 1.13  -1.53   2.17 1.00     2164     1930
##    y1[35]   2.52   2.55 1.16 1.15   0.64   4.37 1.00     1497     1911
##    y1[36]  -1.18  -1.20 1.16 1.18  -3.03   0.76 1.00     1988     1894
##    y1[37]  -0.92  -0.90 1.11 1.14  -2.72   0.91 1.00     1490     1530
##    y1[38]  -0.30  -0.30 1.15 1.15  -2.12   1.63 1.00     1445     1676
##    y1[39]  -2.74  -2.73 1.14 1.10  -4.70  -0.85 1.00     1747     1854
##    y1[40]   0.59   0.57 1.12 1.09  -1.21   2.40 1.00     1892     2012
##    y1[41]   0.27   0.30 1.14 1.08  -1.65   2.12 1.00     1637     1743
##    y1[42]  -0.76  -0.79 1.11 1.13  -2.57   1.07 1.00     2172     2104
##    y1[43]   2.94   2.94 1.11 1.09   1.12   4.78 1.00     2002     1893
##    y1[44]   1.51   1.53 1.10 1.08  -0.34   3.35 1.00     2201     2060
##    y1[45]   2.72   2.73 1.09 1.09   0.93   4.55 1.00     1819     2050
##    y1[46]  -1.39  -1.39 1.11 1.09  -3.25   0.46 1.00     2211     1845
##    y1[47]  -0.62  -0.62 1.11 1.15  -2.42   1.18 1.00     2023     1767
##    y1[48]   2.09   2.10 1.13 1.11   0.24   3.94 1.00     1961     2098
##    y1[49]   1.87   1.86 1.10 1.07   0.08   3.66 1.00     1938     1919
##    y1[50]   2.72   2.72 1.17 1.16   0.81   4.65 1.00     2026     1901
##    m1[1]    0.17   0.17 0.34 0.34  -0.40   0.74 1.00      764     1022
##    m1[2]   -0.52  -0.51 0.29 0.28  -1.02  -0.04 1.00      627      929
##    m1[3]    1.26   1.26 0.31 0.30   0.74   1.77 1.00      818     1072
##    m1[4]   -0.08  -0.07 0.30 0.29  -0.60   0.41 1.00      582      867
##    m1[5]   -0.22  -0.22 0.33 0.31  -0.76   0.32 1.00     2036     1508
##    m1[6]   -0.76  -0.76 0.29 0.28  -1.27  -0.29 1.00      666      979
##    m1[7]    0.34   0.34 0.30 0.29  -0.16   0.82 1.00     1888     1591
##    m1[8]    3.14   3.13 0.38 0.37   2.52   3.78 1.00     1218     1330
##    m1[9]   -0.97  -0.97 0.41 0.40  -1.63  -0.28 1.00      836     1029
##    m1[10]   0.24   0.23 0.39 0.38  -0.38   0.90 1.00      872     1056
##    m1[11]  -0.12  -0.11 0.30 0.29  -0.63   0.38 1.00      584      868
##    m1[12]   1.25   1.26 0.29 0.28   0.77   1.72 1.00     1676     1404
##    m1[13]   2.15   2.16 0.31 0.29   1.64   2.64 1.00     1529     1534
##    m1[14]   1.61   1.60 0.35 0.34   1.04   2.20 1.00      895     1191
##    m1[15]   2.68   2.68 0.35 0.35   2.09   3.24 1.00     1152     1630
##    m1[16]   3.08   3.09 0.36 0.34   2.49   3.67 1.00     1457     1490
##    m1[17]   2.78   2.77 0.37 0.36   2.20   3.40 1.00     1097     1247
##    m1[18]   1.46   1.45 0.36 0.34   0.90   2.05 1.00      881     1226
##    m1[19]  -0.25  -0.26 0.33 0.31  -0.80   0.29 1.00     2049     1453
##    m1[20]   3.57   3.56 0.39 0.38   2.91   4.20 1.00     1454     1394
##    m1[21]  -2.11  -2.10 0.34 0.33  -2.67  -1.54 1.00     1135     1400
##    m1[22]   1.09   1.09 0.31 0.30   0.56   1.61 1.00      794     1068
##    m1[23]   1.25   1.27 0.36 0.35   0.65   1.83 1.00      598      868
##    m1[24]   3.36   3.36 0.38 0.39   2.72   3.99 1.00     1418     1610
##    m1[25]   1.96   1.96 0.32 0.32   1.42   2.48 1.00      941     1319
##    m1[26]  -0.49  -0.49 0.37 0.37  -1.11   0.13 1.00      797      956
##    m1[27]  -0.70  -0.69 0.29 0.28  -1.21  -0.22 1.00      656      987
##    m1[28]   2.18   2.19 0.31 0.29   1.67   2.67 1.00     1526     1534
##    m1[29]   2.36   2.37 0.32 0.30   1.85   2.87 1.00     1511     1484
##    m1[30]   2.53   2.53 0.34 0.35   1.95   3.08 1.00     1104     1592
##    m1[31]  -1.73  -1.72 0.32 0.30  -2.26  -1.21 1.00      959     1187
##    m1[32]  -1.56  -1.57 0.52 0.51  -2.38  -0.68 1.00      998     1302
##    m1[33]   1.78   1.79 0.39 0.40   1.12   2.43 1.00      637      943
##    m1[34]   0.34   0.34 0.30 0.29  -0.15   0.83 1.00     1888     1591
##    m1[35]   2.50   2.49 0.36 0.35   1.93   3.11 1.00     1025     1296
##    m1[36]  -1.18  -1.19 0.39 0.38  -1.83  -0.54 1.00     2236     1548
##    m1[37]  -0.88  -0.88 0.29 0.28  -1.38  -0.41 1.00      690      964
##    m1[38]  -0.29  -0.28 0.29 0.28  -0.80   0.19 1.00      599      862
##    m1[39]  -2.73  -2.72 0.37 0.37  -3.35  -2.10 1.00     1431     1464
##    m1[40]   0.57   0.57 0.33 0.32   0.02   1.11 1.00      766     1101
##    m1[41]   0.28   0.30 0.31 0.30  -0.24   0.78 1.00      568      814
##    m1[42]  -0.77  -0.78 0.36 0.34  -1.37  -0.17 1.00     2187     1456
##    m1[43]   2.94   2.93 0.37 0.36   2.34   3.57 1.00     1165     1276
##    m1[44]   1.51   1.53 0.29 0.27   1.03   1.99 1.00     1618     1569
##    m1[45]   2.70   2.70 0.35 0.35   2.11   3.26 1.00     1158     1630
##    m1[46]  -1.36  -1.36 0.30 0.29  -1.88  -0.87 1.00      822     1052
##    m1[47]  -0.58  -0.58 0.29 0.28  -1.09  -0.11 1.00      637      938
##    m1[48]   2.09   2.09 0.31 0.29   1.58   2.57 1.00     1540     1534
##    m1[49]   1.88   1.89 0.30 0.28   1.38   2.36 1.00     1568     1612
##    m1[50]   2.75   2.74 0.37 0.36   2.17   3.37 1.00     1088     1280

tb=tibble(xa=runif(n,-2,2),xb=runif(n,-2,2),
          ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
          y=rnorm(n,xa+xb-xa*xb+(ca=='b')*2+(cb=='b')*2-(ca=='b')*(cb=='b'),1))

grid.arrange(qplot(data=tb,xa,y,col=ca),
             qplot(data=tb,xa,y,col=cb),
             qplot(data=tb,xb,y,col=ca),
             qplot(data=tb,xb,y,col=cb))

fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
mle

mcmc=goMCMC(mdl,data)

seeMCMC(mcmc,exc='m',ch=F)

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

grid.arrange(
  qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
grid.arrange(
  qplot(data=tb,1:n,y,col=ca)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  qplot(data=tb,1:n,y,col=cb)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  ncol=2
)
}


f0=formula(y~xa+xb+ca+cb)
f1=formula(y~xa+xb+ca*cb)
f2=formula(y~xa*xb+ca*cb)

fn(f0)
## Initial log joint probability = -7137.98 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       22      -43.6763    0.00019513    0.00152004           1           1       36    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -46.54 -46.18 1.88 1.68 -50.15 -44.14 1.00      701      926
##    b[1]     0.42   0.41 0.38 0.38  -0.19   1.04 1.00      671     1000
##    b[2]     0.64   0.64 0.19 0.19   0.32   0.95 1.00     2876     1409
##    b[3]     0.70   0.70 0.21 0.19   0.36   1.04 1.00     2291     1320
##    b[4]     0.84   0.83 0.45 0.46   0.11   1.58 1.00      772      965
##    b[5]     1.94   1.95 0.45 0.43   1.20   2.70 1.01      725     1008
##    s        1.58   1.56 0.18 0.17   1.31   1.90 1.00     2042     1585
##    y1[1]    2.72   2.71 1.72 1.70   0.00   5.56 1.00     1977     1844
##    y1[2]    1.31   1.30 1.63 1.64  -1.34   3.95 1.00     1829     1821
##    y1[3]    1.97   1.99 1.71 1.64  -0.90   4.80 1.00     1681     1965
##    y1[4]    5.10   5.12 1.72 1.65   2.20   7.84 1.00     1928     1801
##    y1[5]    0.15   0.12 1.63 1.60  -2.47   2.85 1.00     1674     1882
##    y1[6]    3.91   3.93 1.64 1.61   1.17   6.57 1.00     2045     1969
##    y1[7]    3.44   3.44 1.64 1.68   0.73   6.10 1.00     2038     1734
##    y1[8]   -0.48  -0.46 1.73 1.67  -3.30   2.30 1.00     1988     1945
##    y1[9]    4.08   4.05 1.70 1.66   1.29   6.96 1.00     2164     1934
##    y1[10]   1.17   1.18 1.66 1.68  -1.53   3.82 1.00     1874     1931
##    y1[11]   1.62   1.61 1.64 1.63  -1.14   4.39 1.00     2033     1995
##    y1[12]   1.49   1.46 1.66 1.63  -1.33   4.30 1.00     1852     1881
##    y1[13]   1.26   1.26 1.64 1.69  -1.43   3.92 1.00     1970     1931
##    y1[14]   2.62   2.63 1.69 1.65  -0.18   5.38 1.00     1995     1857
##    y1[15]   2.22   2.22 1.70 1.67  -0.52   4.93 1.00     1941     1770
##    y1[16]  -0.93  -0.97 1.64 1.60  -3.61   1.83 1.00     1858     1940
##    y1[17]   0.42   0.45 1.65 1.63  -2.31   3.11 1.00     1953     1896
##    y1[18]  -1.00  -0.96 1.69 1.63  -3.84   1.66 1.00     2093     1915
##    y1[19]   2.30   2.33 1.65 1.65  -0.45   4.96 1.00     1895     1997
##    y1[20]   2.40   2.40 1.72 1.65  -0.42   5.24 1.00     2046     1958
##    y1[21]   2.83   2.80 1.68 1.62  -0.05   5.57 1.00     1972     1969
##    y1[22]   2.73   2.73 1.64 1.65   0.08   5.43 1.00     1957     1919
##    y1[23]  -0.87  -0.87 1.66 1.68  -3.54   1.84 1.00     1781     1756
##    y1[24]   1.69   1.71 1.73 1.70  -1.14   4.62 1.00     1878     1496
##    y1[25]   3.87   3.84 1.63 1.60   1.29   6.64 1.00     1822     1972
##    y1[26]   2.21   2.17 1.63 1.60  -0.48   4.90 1.00     1814     1742
##    y1[27]   1.17   1.16 1.60 1.61  -1.41   3.86 1.00     1980     1879
##    y1[28]   3.89   3.88 1.67 1.62   1.10   6.63 1.00     2045     2013
##    y1[29]   0.96   0.92 1.71 1.70  -1.82   3.66 1.00     1835     2012
##    y1[30]   3.16   3.11 1.61 1.60   0.51   5.88 1.00     2156     2058
##    y1[31]   0.61   0.60 1.66 1.60  -2.05   3.34 1.00     2026     1722
##    y1[32]   0.90   0.93 1.67 1.65  -1.84   3.62 1.00     1924     1765
##    y1[33]  -0.46  -0.45 1.65 1.68  -3.16   2.20 1.00     1975     1693
##    y1[34]  -0.96  -1.00 1.69 1.68  -3.67   1.84 1.00     1825     1835
##    y1[35]   3.57   3.59 1.66 1.64   0.88   6.26 1.00     2024     2014
##    y1[36]   2.94   2.94 1.68 1.66   0.13   5.73 1.00     2059     1932
##    y1[37]   0.59   0.62 1.61 1.57  -2.11   3.22 1.00     1897     1857
##    y1[38]   2.02   2.03 1.64 1.53  -0.64   4.74 1.00     1847     1721
##    y1[39]   0.66   0.68 1.68 1.67  -2.11   3.36 1.00     1595     1888
##    y1[40]   2.03   2.03 1.66 1.64  -0.66   4.56 1.00     1898     1916
##    y1[41]   0.14   0.16 1.58 1.61  -2.52   2.62 1.00     1949     1778
##    y1[42]   2.12   2.08 1.66 1.68  -0.53   4.80 1.00     1910     1972
##    y1[43]   3.32   3.35 1.66 1.57   0.62   6.13 1.00     2055     1714
##    y1[44]   3.31   3.35 1.66 1.62   0.59   6.01 1.00     1929     1511
##    y1[45]   4.81   4.81 1.67 1.63   2.09   7.51 1.00     2110     1931
##    y1[46]   3.90   3.88 1.71 1.64   1.03   6.70 1.00     1860     1749
##    y1[47]   2.19   2.24 1.66 1.62  -0.61   4.85 1.00     1660     1902
##    y1[48]   1.52   1.54 1.65 1.59  -1.25   4.19 1.00     1973     1786
##    y1[49]   1.67   1.68 1.67 1.66  -1.13   4.45 1.00     1802     1894
##    y1[50]   2.85   2.88 1.70 1.73   0.05   5.59 1.00     2096     1973
##    m1[1]    2.74   2.74 0.61 0.58   1.72   3.74 1.00     1813     1659
##    m1[2]    1.33   1.35 0.52 0.51   0.46   2.19 1.00     1493     1162
##    m1[3]    1.95   1.95 0.47 0.49   1.20   2.71 1.00      970     1269
##    m1[4]    5.09   5.08 0.54 0.56   4.19   5.96 1.00     1645     1280
##    m1[5]    0.15   0.15 0.38 0.38  -0.47   0.78 1.00      694     1101
##    m1[6]    3.94   3.93 0.52 0.51   3.07   4.81 1.00     1766     1649
##    m1[7]    3.44   3.45 0.58 0.57   2.49   4.37 1.00     1780     1519
##    m1[8]   -0.42  -0.43 0.52 0.51  -1.28   0.47 1.00     1651     1313
##    m1[9]    4.05   4.05 0.51 0.51   3.22   4.90 1.00     1703     1600
##    m1[10]   1.25   1.26 0.45 0.45   0.47   1.97 1.00     1332     1430
##    m1[11]   1.70   1.70 0.47 0.46   0.91   2.43 1.00     1387     1470
##    m1[12]   1.50   1.49 0.43 0.43   0.80   2.21 1.00      816     1178
##    m1[13]   1.28   1.28 0.45 0.45   0.57   2.03 1.00     1082     1205
##    m1[14]   2.64   2.64 0.48 0.46   1.83   3.44 1.00     1495     1622
##    m1[15]   2.28   2.28 0.53 0.50   1.40   3.16 1.00     1896     1337
##    m1[16]  -0.95  -0.95 0.50 0.50  -1.79  -0.14 1.00     1075     1351
##    m1[17]   0.41   0.41 0.43 0.42  -0.30   1.14 1.00     1295     1616
##    m1[18]  -0.99  -0.99 0.49 0.48  -1.79  -0.21 1.00      965      933
##    m1[19]   2.27   2.27 0.46 0.44   1.55   3.03 1.00     1314     1297
##    m1[20]   2.39   2.40 0.70 0.70   1.21   3.51 1.00     1559     1159
##    m1[21]   2.87   2.87 0.51 0.48   2.03   3.71 1.00     1824     1265
##    m1[22]   2.76   2.78 0.49 0.49   1.95   3.55 1.00     1479     1285
##    m1[23]  -0.91  -0.91 0.48 0.49  -1.69  -0.13 1.00      985     1257
##    m1[24]   1.67   1.67 0.61 0.57   0.70   2.72 1.00     1433     1169
##    m1[25]   3.83   3.83 0.49 0.48   3.06   4.62 1.00     1427     1680
##    m1[26]   2.23   2.23 0.50 0.52   1.42   3.04 1.00     1127     1283
##    m1[27]   1.17   1.18 0.42 0.41   0.46   1.88 1.00     1268     1124
##    m1[28]   3.90   3.91 0.51 0.49   3.04   4.73 1.00     1715     1127
##    m1[29]   0.96   0.96 0.57 0.56   0.02   1.91 1.01      924      922
##    m1[30]   3.20   3.20 0.49 0.45   2.40   4.00 1.00     1733     1130
##    m1[31]   0.53   0.54 0.48 0.44  -0.25   1.32 1.00     1389     1592
##    m1[32]   0.88   0.89 0.44 0.42   0.15   1.60 1.00     1288     1595
##    m1[33]  -0.50  -0.50 0.45 0.45  -1.25   0.26 1.00      969     1414
##    m1[34]  -0.98  -0.99 0.57 0.56  -1.91  -0.06 1.00     1028     1219
##    m1[35]   3.60   3.60 0.65 0.64   2.57   4.67 1.00     1387     1468
##    m1[36]   2.92   2.94 0.45 0.42   2.18   3.65 1.00     1574     1131
##    m1[37]   0.54   0.53 0.37 0.38  -0.06   1.16 1.00      671     1006
##    m1[38]   2.04   2.03 0.42 0.43   1.37   2.73 1.00     1150     1308
##    m1[39]   0.69   0.68 0.60 0.58  -0.27   1.70 1.00     1596     1253
##    m1[40]   2.01   2.00 0.51 0.49   1.18   2.85 1.00     1783     1219
##    m1[41]   0.06   0.07 0.41 0.43  -0.60   0.75 1.00      879     1292
##    m1[42]   2.08   2.08 0.43 0.43   1.39   2.77 1.00     1194     1309
##    m1[43]   3.39   3.41 0.58 0.58   2.45   4.34 1.00     1800     1478
##    m1[44]   3.33   3.34 0.45 0.41   2.60   4.06 1.00     1545     1103
##    m1[45]   4.83   4.82 0.61 0.61   3.83   5.82 1.00     1916     1609
##    m1[46]   3.89   3.88 0.52 0.51   3.04   4.73 1.00     1255     1424
##    m1[47]   2.18   2.18 0.41 0.42   1.51   2.86 1.00     1121     1325
##    m1[48]   1.48   1.49 0.46 0.45   0.71   2.24 1.00     1348      952
##    m1[49]   1.64   1.64 0.54 0.54   0.76   2.52 1.01      926      886
##    m1[50]   2.79   2.80 0.54 0.53   1.90   3.68 1.00     1650     1426

fn(f1)
## Initial log joint probability = -3972.6 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       27      -43.0984   0.000110501   0.000793672           1           1       40    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -46.49 -46.14 1.94 1.77 -50.18 -44.06 1.00      832     1244
##    b[1]     0.26   0.25 0.42 0.42  -0.45   0.97 1.00      856     1184
##    b[2]     0.58   0.59 0.19 0.20   0.26   0.89 1.00     2297     1670
##    b[3]     0.70   0.70 0.20 0.19   0.39   1.03 1.00     2955     1718
##    b[4]     1.22   1.22 0.60 0.58   0.26   2.20 1.01      705      910
##    b[5]     2.39   2.38 0.67 0.68   1.30   3.47 1.01      596      790
##    b[6]    -0.94  -0.96 0.91 0.94  -2.41   0.54 1.00      501      689
##    s        1.58   1.56 0.18 0.17   1.33   1.89 1.00     2207     1792
##    y1[1]    2.88   2.90 1.71 1.69   0.10   5.72 1.00     2016     2010
##    y1[2]    1.68   1.70 1.64 1.62  -0.91   4.34 1.00     1802     1852
##    y1[3]    1.71   1.73 1.69 1.64  -1.04   4.34 1.00     1888     1895
##    y1[4]    4.77   4.72 1.66 1.70   2.09   7.55 1.00     1536     1727
##    y1[5]   -0.01  -0.04 1.63 1.57  -2.67   2.67 1.00     2055     1877
##    y1[6]    4.10   4.13 1.70 1.65   1.32   6.94 1.00     2040     1931
##    y1[7]    3.63   3.65 1.71 1.70   0.77   6.42 1.00     1861     1865
##    y1[8]   -0.17  -0.14 1.65 1.59  -2.91   2.49 1.00     1795     1701
##    y1[9]    4.21   4.16 1.68 1.62   1.46   6.93 1.00     1787     1687
##    y1[10]   1.49   1.45 1.64 1.62  -1.23   4.16 1.00     1954     1826
##    y1[11]   1.83   1.81 1.66 1.67  -0.87   4.52 1.00     1956     1856
##    y1[12]   1.29   1.29 1.68 1.66  -1.52   4.04 1.00     1895     2102
##    y1[13]   1.04   1.07 1.73 1.71  -1.81   3.82 1.00     1668     1575
##    y1[14]   2.86   2.82 1.67 1.64   0.18   5.68 1.00     1729     1832
##    y1[15]   2.11   2.08 1.66 1.66  -0.63   4.77 1.00     2175     1836
##    y1[16]  -1.09  -1.07 1.69 1.68  -3.89   1.62 1.00     1758     1403
##    y1[17]   0.67   0.63 1.64 1.68  -1.89   3.33 1.00     1927     1784
##    y1[18]  -1.13  -1.12 1.65 1.63  -3.92   1.51 1.00     2003     2009
##    y1[19]   2.03   2.06 1.70 1.70  -0.75   4.74 1.00     1987     1917
##    y1[20]   2.77   2.73 1.76 1.66  -0.01   5.65 1.00     1970     1915
##    y1[21]   2.74   2.68 1.65 1.57   0.16   5.45 1.00     1969     1826
##    y1[22]   2.91   2.90 1.71 1.68   0.18   5.64 1.00     1939     1729
##    y1[23]  -1.04  -1.05 1.64 1.60  -3.74   1.63 1.00     1657     1820
##    y1[24]   1.39   1.42 1.73 1.73  -1.50   4.20 1.00     2005     1954
##    y1[25]   4.07   4.05 1.67 1.69   1.32   6.76 1.00     1722     1577
##    y1[26]   1.94   1.99 1.69 1.64  -0.95   4.78 1.00     1982     1704
##    y1[27]   1.42   1.39 1.68 1.63  -1.35   4.33 1.00     2016     1974
##    y1[28]   3.69   3.70 1.70 1.62   0.85   6.54 1.00     2003     1706
##    y1[29]   0.81   0.82 1.67 1.63  -1.94   3.60 1.00     1836     1527
##    y1[30]   2.99   3.00 1.68 1.69   0.26   5.78 1.00     1840     1852
##    y1[31]   0.72   0.71 1.63 1.59  -1.95   3.41 1.00     1949     1850
##    y1[32]   1.06   1.10 1.65 1.62  -1.73   3.73 1.00     1661     1908
##    y1[33]  -0.67  -0.68 1.63 1.60  -3.30   2.00 1.00     1981     1547
##    y1[34]  -1.05  -1.07 1.71 1.69  -3.87   1.83 1.00     1969     2001
##    y1[35]   3.17   3.14 1.76 1.73   0.25   6.05 1.00     1923     1782
##    y1[36]   2.72   2.71 1.68 1.65   0.00   5.48 1.00     1980     1684
##    y1[37]   0.39   0.40 1.66 1.65  -2.39   3.09 1.00     1924     1645
##    y1[38]   2.34   2.34 1.69 1.75  -0.38   5.08 1.00     1747     1669
##    y1[39]   1.03   1.04 1.69 1.67  -1.69   3.87 1.00     1836     1855
##    y1[40]   1.91   1.88 1.68 1.71  -0.85   4.67 1.00     2131     1973
##    y1[41]  -0.11  -0.11 1.62 1.65  -2.64   2.58 1.00     1933     2049
##    y1[42]   2.36   2.31 1.66 1.62  -0.33   5.13 1.00     2027     1976
##    y1[43]   3.51   3.53 1.69 1.67   0.65   6.21 1.00     2131     2058
##    y1[44]   3.11   3.09 1.67 1.68   0.31   5.82 1.00     1901     1894
##    y1[45]   5.05   5.03 1.74 1.69   2.19   7.96 1.00     1921     1776
##    y1[46]   3.50   3.50 1.71 1.69   0.67   6.23 1.00     1831     1895
##    y1[47]   2.43   2.39 1.70 1.60  -0.38   5.31 1.00     1818     1893
##    y1[48]   1.74   1.75 1.65 1.59  -1.04   4.42 1.00     1942     1933
##    y1[49]   1.53   1.52 1.63 1.62  -1.12   4.16 1.00     1924     1787
##    y1[50]   2.86   2.87 1.69 1.64   0.11   5.66 1.00     1753     2015
##    m1[1]    2.91   2.92 0.62 0.58   1.87   3.95 1.00     2064     1603
##    m1[2]    1.63   1.63 0.58 0.60   0.70   2.58 1.00     1311     1475
##    m1[3]    1.74   1.74 0.52 0.53   0.86   2.60 1.00      860     1350
##    m1[4]    4.79   4.79 0.59 0.59   3.80   5.76 1.00     1743     1593
##    m1[5]    0.00   0.00 0.42 0.42  -0.70   0.71 1.00      906     1103
##    m1[6]    4.13   4.13 0.56 0.53   3.23   5.04 1.00     1670     1740
##    m1[7]    3.57   3.58 0.59 0.61   2.58   4.53 1.00     1858     1304
##    m1[8]   -0.13  -0.13 0.60 0.62  -1.12   0.86 1.00     1243     1447
##    m1[9]    4.26   4.25 0.56 0.55   3.38   5.16 1.00     1516     1591
##    m1[10]   1.42   1.42 0.49 0.50   0.63   2.22 1.00     1294     1480
##    m1[11]   1.85   1.84 0.49 0.49   1.03   2.67 1.00     1383     1398
##    m1[12]   1.31   1.31 0.47 0.48   0.53   2.09 1.00      838     1270
##    m1[13]   1.04   1.04 0.51 0.51   0.17   1.88 1.00      716     1095
##    m1[14]   2.86   2.86 0.53 0.52   1.99   3.74 1.00     1309     1229
##    m1[15]   2.12   2.11 0.54 0.54   1.23   3.00 1.00     2431     1565
##    m1[16]  -1.10  -1.10 0.53 0.51  -1.97  -0.21 1.00     1281     1477
##    m1[17]   0.66   0.66 0.50 0.52  -0.15   1.46 1.00     1128     1311
##    m1[18]  -1.11  -1.11 0.51 0.52  -1.94  -0.27 1.00     1396     1722
##    m1[19]   2.02   2.03 0.51 0.51   1.17   2.84 1.00     1746     1426
##    m1[20]   2.81   2.77 0.80 0.81   1.55   4.14 1.00     1097     1248
##    m1[21]   2.69   2.68 0.52 0.50   1.84   3.55 1.00     2327     1655
##    m1[22]   2.95   2.94 0.52 0.52   2.10   3.81 1.00     1516     1397
##    m1[23]  -1.05  -1.05 0.51 0.50  -1.86  -0.20 1.00     1294     1575
##    m1[24]   1.40   1.41 0.65 0.65   0.29   2.44 1.00     1888     1492
##    m1[25]   4.09   4.07 0.56 0.56   3.20   5.01 1.00     1221     1464
##    m1[26]   2.00   2.01 0.56 0.57   1.03   2.92 1.00      856     1405
##    m1[27]   1.44   1.43 0.49 0.51   0.64   2.25 1.00     1133     1533
##    m1[28]   3.68   3.68 0.53 0.52   2.82   4.54 1.00     2077     1739
##    m1[29]   0.86   0.85 0.56 0.57  -0.07   1.80 1.00     1453     1356
##    m1[30]   3.00   3.00 0.50 0.49   2.15   3.84 1.00     2153     1742
##    m1[31]   0.73   0.72 0.52 0.54  -0.11   1.57 1.00     1296     1433
##    m1[32]   1.07   1.06 0.48 0.49   0.29   1.87 1.00     1208     1514
##    m1[33]  -0.67  -0.68 0.49 0.47  -1.48   0.18 1.00     1076     1188
##    m1[34]  -1.03  -1.03 0.58 0.60  -1.97  -0.08 1.00     1731     1731
##    m1[35]   3.21   3.19 0.72 0.72   2.04   4.38 1.00     1503     1584
##    m1[36]   2.72   2.72 0.48 0.47   1.92   3.52 1.00     1988     1651
##    m1[37]   0.38   0.37 0.42 0.42  -0.32   1.08 1.00      833     1175
##    m1[38]   2.33   2.33 0.52 0.53   1.47   3.22 1.01      919     1267
##    m1[39]   0.99   0.98 0.67 0.64  -0.05   2.13 1.00     1162     1407
##    m1[40]   1.83   1.83 0.53 0.53   0.94   2.70 1.00     2273     1568
##    m1[41]  -0.13  -0.13 0.46 0.45  -0.89   0.67 1.00      891     1034
##    m1[42]   2.36   2.36 0.52 0.52   1.52   3.23 1.01      961     1213
##    m1[43]   3.51   3.53 0.59 0.61   2.52   4.47 1.00     1841     1422
##    m1[44]   3.11   3.12 0.48 0.46   2.32   3.93 1.00     1915     1623
##    m1[45]   5.04   5.02 0.65 0.65   4.02   6.07 1.00     1727     1725
##    m1[46]   3.54   3.53 0.60 0.60   2.53   4.52 1.00     1393     1529
##    m1[47]   2.48   2.47 0.52 0.53   1.64   3.36 1.01      903     1239
##    m1[48]   1.75   1.75 0.52 0.53   0.92   2.60 1.00     1224     1395
##    m1[49]   1.50   1.49 0.55 0.56   0.58   2.38 1.00     1202     1354
##    m1[50]   2.91   2.92 0.55 0.55   1.99   3.82 1.00     1617     1334

fn(f2)
## Initial log joint probability = -126.472 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       31      -24.9459   4.85521e-05   0.000907066      0.6786      0.6786       35    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -29.35 -29.00 2.14 2.01 -33.42 -26.54 1.00      792     1205
##    b[1]     0.51   0.51 0.30 0.28   0.03   1.00 1.01      709     1107
##    b[2]     0.92   0.92 0.15 0.15   0.67   1.17 1.00     1834     1212
##    b[3]     0.79   0.79 0.14 0.14   0.56   1.03 1.00     3182     1847
##    b[4]     1.15   1.15 0.43 0.43   0.46   1.88 1.01      511      843
##    b[5]     1.98   1.99 0.46 0.45   1.24   2.73 1.01      452      794
##    b[6]    -0.89  -0.88 0.14 0.14  -1.12  -0.65 1.00     2410     1681
##    b[7]    -1.04  -1.01 0.66 0.66  -2.15  -0.02 1.01      394      596
##    s        1.11   1.10 0.13 0.12   0.93   1.34 1.00     2082     1563
##    y1[1]    5.50   5.48 1.27 1.25   3.51   7.61 1.00     1852     2011
##    y1[2]    2.99   3.01 1.23 1.21   0.95   5.03 1.00     2075     1933
##    y1[3]    1.33   1.35 1.18 1.14  -0.65   3.25 1.00     1959     1648
##    y1[4]    3.50   3.52 1.18 1.16   1.53   5.35 1.00     2097     1857
##    y1[5]    0.17   0.18 1.14 1.08  -1.68   2.08 1.00     1678     1735
##    y1[6]    3.76   3.76 1.23 1.17   1.66   5.86 1.00     1931     1933
##    y1[7]    2.07   2.11 1.21 1.20   0.05   4.03 1.00     2000     2059
##    y1[8]   -1.90  -1.88 1.22 1.19  -3.94   0.14 1.00     2094     1892
##    y1[9]    3.45   3.44 1.17 1.16   1.55   5.40 1.00     2050     2012
##    y1[10]   2.27   2.27 1.20 1.23   0.31   4.24 1.00     1767     1917
##    y1[11]   2.84   2.87 1.16 1.11   0.94   4.79 1.00     1900     1906
##    y1[12]   1.34   1.32 1.18 1.15  -0.59   3.33 1.00     1772     1934
##    y1[13]   2.02   2.01 1.19 1.21   0.11   3.95 1.00     1905     1968
##    y1[14]   3.74   3.75 1.19 1.18   1.83   5.64 1.00     1950     2014
##    y1[15]   1.93   1.92 1.19 1.21   0.00   3.87 1.00     1838     1821
##    y1[16]  -1.65  -1.66 1.16 1.12  -3.53   0.28 1.00     1740     2148
##    y1[17]   0.21   0.18 1.14 1.09  -1.72   2.13 1.00     2043     1940
##    y1[18]  -2.29  -2.27 1.20 1.15  -4.27  -0.27 1.00     1732     2010
##    y1[19]   1.00   1.02 1.19 1.20  -0.91   2.97 1.00     1956     2082
##    y1[20]   5.28   5.30 1.31 1.27   3.16   7.47 1.00     1979     1831
##    y1[21]   3.28   3.29 1.17 1.14   1.30   5.16 1.00     2040     1918
##    y1[22]   2.58   2.62 1.17 1.16   0.66   4.50 1.00     1642     1816
##    y1[23]  -1.84  -1.84 1.16 1.15  -3.68   0.06 1.00     1878     1947
##    y1[24]   0.38   0.39 1.18 1.15  -1.56   2.29 1.00     1816     1973
##    y1[25]   3.38   3.40 1.17 1.14   1.36   5.28 1.00     1855     2059
##    y1[26]   1.26   1.27 1.19 1.11  -0.71   3.23 1.00     1747     1912
##    y1[27]   1.80   1.83 1.19 1.19  -0.24   3.68 1.00     2050     2052
##    y1[28]   4.39   4.43 1.21 1.20   2.41   6.33 1.00     1863     1931
##    y1[29]   2.24   2.25 1.19 1.17   0.34   4.18 1.00     1813     1844
##    y1[30]   3.59   3.60 1.19 1.19   1.67   5.61 1.00     2083     2056
##    y1[31]   1.04   1.03 1.17 1.18  -0.88   2.93 1.00     1968     1888
##    y1[32]   1.45   1.46 1.19 1.22  -0.52   3.34 1.00     1777     1932
##    y1[33]  -0.51  -0.50 1.14 1.08  -2.36   1.36 1.00     1815     1972
##    y1[34]  -1.77  -1.78 1.19 1.15  -3.76   0.16 1.00     1862     1941
##    y1[35]   5.45   5.46 1.28 1.28   3.40   7.53 1.00     1951     1892
##    y1[36]   2.70   2.69 1.16 1.13   0.83   4.67 1.00     1801     1933
##    y1[37]   0.65   0.64 1.17 1.15  -1.27   2.56 1.00     1594     1682
##    y1[38]   2.13   2.14 1.22 1.19   0.11   4.11 1.00     1953     2088
##    y1[39]  -0.33  -0.37 1.20 1.17  -2.32   1.65 1.00     1788     1887
##    y1[40]   0.71   0.72 1.21 1.19  -1.26   2.69 1.00     1854     1896
##    y1[41]   0.56   0.54 1.19 1.18  -1.41   2.52 1.00     1950     1950
##    y1[42]   2.29   2.29 1.18 1.14   0.33   4.25 1.00     1918     2009
##    y1[43]   2.25   2.21 1.20 1.15   0.28   4.28 1.00     1832     1892
##    y1[44]   3.28   3.27 1.18 1.18   1.34   5.21 1.00     1984     1932
##    y1[45]   2.63   2.66 1.30 1.34   0.49   4.69 1.00     1897     1748
##    y1[46]   4.09   4.08 1.21 1.16   2.09   6.12 1.00     1670     1831
##    y1[47]   2.21   2.21 1.19 1.16   0.28   4.20 1.00     1812     1951
##    y1[48]   2.50   2.48 1.20 1.17   0.54   4.48 1.00     2010     1973
##    y1[49]   1.96   1.94 1.20 1.17   0.02   4.01 1.00     1899     1916
##    y1[50]   2.98   2.99 1.17 1.18   1.03   4.85 1.00     1940     1882
##    m1[1]    5.53   5.52 0.60 0.60   4.56   6.51 1.00     2471     1644
##    m1[2]    2.97   2.96 0.47 0.45   2.19   3.75 1.00     1371     1617
##    m1[3]    1.31   1.30 0.37 0.37   0.69   1.91 1.01      768     1431
##    m1[4]    3.52   3.51 0.49 0.48   2.74   4.32 1.00     1684     1726
##    m1[5]    0.15   0.15 0.30 0.28  -0.34   0.65 1.01      755      963
##    m1[6]    3.76   3.76 0.39 0.37   3.09   4.40 1.00     1536     1573
##    m1[7]    2.13   2.13 0.47 0.47   1.33   2.90 1.00     1843     1649
##    m1[8]   -1.92  -1.92 0.52 0.52  -2.78  -1.06 1.00     1448     1301
##    m1[9]    3.46   3.46 0.41 0.40   2.77   4.12 1.00     1578     1654
##    m1[10]   2.28   2.28 0.37 0.37   1.68   2.89 1.00     1439     1143
##    m1[11]   2.84   2.84 0.38 0.39   2.22   3.47 1.00     1525     1416
##    m1[12]   1.32   1.32 0.33 0.33   0.77   1.87 1.01      720     1193
##    m1[13]   2.06   2.05 0.39 0.38   1.42   2.71 1.00      724     1347
##    m1[14]   3.81   3.81 0.40 0.39   3.11   4.44 1.00     1398     1654
##    m1[15]   1.94   1.94 0.38 0.36   1.33   2.55 1.00     2245     1679
##    m1[16]  -1.65  -1.65 0.38 0.38  -2.28  -1.02 1.00     1276     1522
##    m1[17]   0.24   0.25 0.37 0.36  -0.37   0.84 1.00     1111     1115
##    m1[18]  -2.25  -2.25 0.40 0.40  -2.93  -1.60 1.00     1495     1565
##    m1[19]   0.99   1.00 0.38 0.39   0.35   1.58 1.00     1717     1575
##    m1[20]   5.29   5.31 0.70 0.71   4.15   6.45 1.00     1281     1264
##    m1[21]   3.30   3.29 0.38 0.37   2.66   3.93 1.00     2237     1543
##    m1[22]   2.55   2.55 0.37 0.36   1.96   3.16 1.00     1555     1590
##    m1[23]  -1.85  -1.85 0.38 0.38  -2.48  -1.23 1.00     1315     1580
##    m1[24]   0.37   0.38 0.47 0.49  -0.42   1.11 1.00     1819     1575
##    m1[25]   3.32   3.32 0.41 0.41   2.64   4.00 1.00     1242     1529
##    m1[26]   1.26   1.26 0.41 0.42   0.57   1.93 1.01      807     1516
##    m1[27]   1.81   1.81 0.36 0.36   1.22   2.42 1.00     1152     1224
##    m1[28]   4.41   4.39 0.40 0.40   3.74   5.07 1.00     2143     1790
##    m1[29]   2.23   2.23 0.46 0.45   1.48   2.96 1.00     1297     1170
##    m1[30]   3.59   3.58 0.37 0.37   2.97   4.20 1.00     2141     1695
##    m1[31]   1.07   1.07 0.38 0.38   0.45   1.70 1.00     1279      955
##    m1[32]   1.49   1.49 0.35 0.36   0.90   2.07 1.00     1207      897
##    m1[33]  -0.49  -0.49 0.35 0.34  -1.05   0.08 1.00      956     1336
##    m1[34]  -1.79  -1.79 0.43 0.44  -2.51  -1.09 1.00     1579     1471
##    m1[35]   5.42   5.42 0.61 0.60   4.44   6.44 1.00     1654     1460
##    m1[36]   2.70   2.70 0.33 0.32   2.14   3.24 1.00     1905     1607
##    m1[37]   0.67   0.66 0.30 0.28   0.18   1.15 1.01      693     1093
##    m1[38]   2.10   2.10 0.37 0.36   1.50   2.69 1.00      844     1261
##    m1[39]  -0.32  -0.32 0.51 0.51  -1.15   0.52 1.00     1188     1482
##    m1[40]   0.74   0.74 0.40 0.40   0.08   1.38 1.00     2184     1654
##    m1[41]   0.59   0.59 0.34 0.33   0.03   1.16 1.00      795     1268
##    m1[42]   2.29   2.30 0.36 0.36   1.69   2.87 1.00      875     1288
##    m1[43]   2.30   2.29 0.46 0.46   1.53   3.04 1.00     1781     1568
##    m1[44]   3.32   3.31 0.34 0.33   2.76   3.87 1.00     1867     1598
##    m1[45]   2.64   2.65 0.60 0.60   1.66   3.61 1.00     1988     1715
##    m1[46]   4.10   4.10 0.43 0.43   3.38   4.80 1.00     1209     1663
##    m1[47]   2.26   2.26 0.36 0.35   1.66   2.86 1.00      827     1254
##    m1[48]   2.50   2.50 0.39 0.39   1.85   3.14 1.00     1262     1437
##    m1[49]   1.95   1.95 0.40 0.39   1.30   2.62 1.00     1040     1273
##    m1[50]   3.00   3.00 0.38 0.39   2.38   3.65 1.00     1577     1571